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INTRODUCTION 



The USAF Western Space and Missile Center is located along the 
southern California coast at Vandenberg AFB, where large Atlas and 
Titan rockets are used to launch polar orbiting satellites. Since 
the rockets are fueled by tanker trucks, USAF has directed much 
effort toward avoiding or mitigating the potential for large, 
accidental, ground level releases of toxic fumes. 

If released, such fumes would disperse downwind according to the 
detailed mechanics of the release, nature of the propellants, and 
prevailing atmospheric conditions and terrain. Even during normal 
launch, massive amounts of propellants are released into the 
atmospheric boundary layer via the rocket exhaust cloud. Thus, 
USAF focuses part of its mitigation effort on field studies and 
computer forecast/nowcasts of such cold and hot spills to prepare 
toxic hazard corridors (THCs) for base and civilian personnel. 

For example, Battelle Corp. conducted the Mt. Iron dispersion study 
from the three south base launch complexes, SLC-4, SLC-5, and 
SLC-6, with 252 zinc sulfide tracer releases from Dec. 1965 to July 
1967 (Hinds and Nickola, 1967). Until the recent Lompoc Valley 
Diffusion Study (LVDE) of potential emissions from the Hypergolic 
Stockpile and Storage Facility (HSSF) (Skupniewicz , et al. 1990, 
1991), Mt. Iron was the only field study involving tracer releases 
conducted from South Vandenberg. Eighty-eighc. releases were 
conducted from SLC-4, 25 from SLC-5, and 139 from SLC-6. 

since SLC-6 is presently decommissioned, SLC-4 is now the main 
launch pad for large Vandenberg rockets. Thus, in an effort to 
assess possible replacements for the currently used Ocean Breeze/ 
Dry Gulch (OBDG) type regression equations (Haugen and Taylor, 
1963), more physically based computer models such as AFTOX and 
HOTMAC/RAPTAD (Kunkel and Izumi, 1991; Yamada and Bunker, 1991) 
have been used to simulate Mt. Iron releases from SLC-4. The Naval 
Postgraduate School (NPS) has made available an eight case subset 
of the Mt. Iron data set to the Vandenberg modeling community and 
has made efforts to coordinate the various simulation studies. 
This report describes the use of LINCOM/RIMPUFF, a dispersion model 
from RISO National Laboratory of Denmark, to simulate the above 
eight representative cases. 



2. MOUNTAIN IRON DATA DESCRIPTION 



2.1 The Vandenberg AFB Domain 

Vandenberg AFB is located about 100 miles northwest of Los Angeles 
on the California coast. The Vandenberg domain features generally 
hilly terrain, punctuated by the east-west oriented Santa Ynez 
River Valley which separates North and South Vandenberg. The South 
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Vandenberg terrain is generally rougher than that of North 
Vandenberg with ridge to canyon elevation differences ranging from 
around 100 to 400m. At Pt. Arguello (near SLC-6) the coastline 
slants sharply eastward away from the generally north-south 
direction. This is repeated further south beyond Vandenberg by a 
similar coastline turning feature at Pt. Concepcion. Honda Ridge 
on South Vandenberg has a high point at Tranquillon Peak of over 
600m and an average height of about 200m. It represents the 
northern edge of the coastal Santa Ynez Mountains. At an average 
elevation of about 60m, Honda Canyon, roughly 2km wide, separates 
Honda Ridge from Target Ridge, the first ridge south of SLC-4 (see 
fig. 1) . 

The climatological flow is from the northwest, controlled by 
subsiding anti-cyclonic flow around the chronic eastern Pacific 
High. Adiabatic compression of this subsiding air mass creates the 
typically strong, low California coastal subsidence inversion seen 
at Vandenberg. The northwesterly flow is augmented by the summer 
time California coast-to-Central-Valley monsoon, as well as the 
local Seabreeze. However, the turning coastline at Pts. Arguello 
and Concepcion induces a local divergence of the climatological 
northwesterlies . Moreover, Honda Ridge, together with the 
typically low subsidence inversion, tends to divert daytime 
northwesterly flows to the west, so that winds over SLC-6 are 
usually channeled to a narrow range about the north-south direction 
(see Vandenberg Meteorology and Plume Dispersion Handbook, Kamada 
et al, 1989, hereafter referred to as the Handbook). 



2.2 General Data Description 

113 tests were conducted at SLC-4 and SLC-5 during Phase I of Mt. 
Iron in the first half of 1966. In the experimental design, 
fluorescent zinc sulfide tracer particles with a mean solid 
diameter of 2.5 nm was released as a wet aerosol fog from a test 
stand at a height of 12 ft, usually for 30 minutes. Sample tracer 
dosages were collected on membrane filters and assayed for 
fluorescence by alpha emission, using a Rankin counter. Due to the 
rugged terrain, the sampling lines were limited to existing paved 
and gravel roads as shown in fig. 2. Sampler spacing ranged between 
approximately 160 - 640 meters (0.1 - 0.4 miles), depending on 
distance from the release point. The 160m spacing was used for 
roads near the release site. A 320m spacing was used along the 
eastern section of Santa Ynez, Honda Canyon, Honda Ridge, and 
Tranquillon Mountain Rds. The 640m spacing was reserved for the 
section of the Sudden Coast Rd . from Pt. Arguello stretching south 
by southeast to Jalama Beach . 

For this data/model comparison, we selected eight representative 
cases. The cases were chosen on the basis of completeness of the 
data set collected by Battelle during Mt. Iron and the significance 
of the presumed flow type as defined in the Handbook. 
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Fig. 1 Terrain Map of South Vandenberg. Inner rectangle shows 
the 100m inner grid area for RIMPUFF. Elevations are in meters. 
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Fig. 2 Sampling routes over South Vandenberg for Mt. Iron. Inner 
rectangle shows the 100m inner grid area for LINCOM. Weather 
towers are numbered. 
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The available Mt. Iron sensor array consisted of towers, sondes, 
and bag samplers. There were: 1) wind speeds, directions, and 
direction standard deviations from 19 towers over 7.5, 30, 60, and 
90 minute averages (shown in Tables 2.1.1, 2.1.2 and fig. 2), 2) 
tethersonde temperatures at approximately 20 meter intervals up to 
260 m. The tethersonde (often referred to as a wiresonde) was 
located at the release site, either VIP-l at the northwest corner 
of SLC-4 or Area 529 at the southwest corner of SLC-5. 3) 
twice-a-day NWS rawinsondes from Building 22, 4) supplemental rawin 
or radiosondes from VIP-1, Scout D, upper Honda Canyon, and the 
Boathouse. The rawinsondes were released as often as every hour 
for 3-4 hours, usually commencing one hour before tracer release. 
5) bag sampling lines were deployed along Mesa, Coast, Surf, 
Arguello, Santa Ynez, Bear Creek, Folo, Tank, Target Ridge, Honda 
Canyon, Tranquillon Mt., Honda Ridge, and Sudden Coast roads, using 
a total of " 330 bag samplers. The Queen Air aircraft, used for 24 
flights during Phase II, was deployed only for case 110 of the 
eight cases selected for study. For case 110 the plume centerlines 
derived from the aircraft data were essentially identical with 
those from ground sampling. 



Table 2.1.1 Tower/Sonde sites 



Program Index. VAFB Site #. 



1 


301 






2 


056 






3 


014 






4 


012 






5 


101 






6 


Target Ridge 




7 


Oil 






8 


300 






9 


054 






10 


100 






11 


VHF 






12 


010 






13 


009 






14 


055 






15 


057 


(Boathouse, BH) 


* 


16 


013 


(Honda Canyon, HC) 


★ 


17 


200 


(Scout-D, SD) 


★ 


18 


VIP-1 




★ 


19 


B22 


(Bldg. 22) 


* 



The asterisks denote rawinsonde launches at the same site as the 
tower 
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Table 2.1.2 Wind Speed, Direction, 0, and 
Standard Deviation, Oq, for Main Sensor Sites 



Site 


28 


31 


48 


VIP-1 




6.7 


1.5 




285 


330 


318 




17 


5 




VHF 


3.8 


7 . 4 






294 


341 






7 


13 




Target 


1 . 2 


9 . 6 


2 . 3 


Ridge 


269 


4 


55 




3 


19 


6 


WTlOl 


4 . 1 


7 . 4 


4 . 1 




300 


330 


325 




8 


15 


8 


Telemetry 


2 . 4 


4 . 1 




Peak 


321 


335 






5 


8 




Ion 


2 . 5 


2.0 


1.0 


Sounder 


259 


328 


235 


Honda Cyn 


5 


4 


2 


LaSalle 


2 . 4 




1 . 0 


Canyon Rd 


328 




212 


Honda Cyn 


5 






WT200 


1.9 


2 . 3 






274 


220 


10 




4 


5 




WT3 01 


2 . 2 








339 


339 






5 


15 




Boat 


3.5 


7 . 5 


7 . 0 


House 


280 


350 


315 




7 


17 


13 



87 


90 


91 


110 


Case 


5 . 2 


3 . 9 


3 . 3 


2.7 


V m/s 


250 


325 


360 


270 


0 deg 


6 


8 


6 


15 


Oq deg 


4 . 4 


5.1 


2 . 8 


5.0 


V 


275 


303 


340 


296 


e 


11 


5 


5 




^0 


1.1 


5 . 8 


6 . 8 


2 . 0 


V 


271 


298 


324 


323 


0 


12 


9 


6 


9 


V 


251 


300 


315 


325 


0 


2 


12 


14 


4 


^0 


2 . 0 


6.0 


5.7 


2.0 


V 


318 


310 


313 


16 


0 


4 


12 


12 


4 


^0 


3 . 4 


2 . 9 


3 . 5 


4 . 1 


V 


256 


1 


360 


249 


0 


7 


6 


7 


8 


^0 


2 . 6 


4 . 0 


2 . 6 


2 . 1 


V 


268 


315 


316 


275 


0 


5 


8 


5 


4 


^0 




1.2 


0.8 


2 . 7 


V 


254 


184 


105 


259 


0 


8 


2 


2 


5 


^0 


3 . 9 


7 . 8 


7 . 6 


2 . 8 


V 


204 


1 


15 


329 


0 


8 


16 


15 


6 


^0 


2 . 3 


13 . 0 


14 . 5 


4 . 6 


V 


213 


309 


313 


239 


0 


7 


26 


29 


9 


^0 



55 

5.0 

321 

7.9 

333 

14 

5.7 

339 

11 

7.4 

303 

15 

4 . 9 

309 

10 

2.9 

288 

6 

2 . 9 

298 

9 

0.5 

256 

7 

7 . 4 

340 

15 

340 
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2.3 Case Study Selection 



As the Mt. Iron study gained momentum, hourly rawinsondes became 
available at up to four locations. However, this was not true in 
the early going, and the wiresonde at the release site was not 
operative for the first two dozen or so cases. Thus, upper air 
data during the early stages of the study was rather sparse. Even 
some of the ground based towers were not yet installed. So, 
although the Mt. Iron study began in the beginning of December 
1965, our first study case with fairly complete data was case 28 
which took place in February of 1966. 

Many of the sensors also were inoperative at various times, 
particularly the wiresonde. Since release site winds and 
temperature structure are usually the most critical inputs when 
simulating plume transport, we tried to limit our selections to 
those cases where wiresonde data was available. 

Within the available data, our other primary consideration was to 
sample as many flow conditions as possible. Table 2.2.1 summarizes 
ten significant flow types, as characterized in the Handbook. 
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Table 2.2.1 Surface Flow Types and Typical Conditions 





FLOW 

TYPE 


FLOW 

ALOFT 

052 


ANGLE 

AT 


INVERSION 

HEIGHT 


USUALLY 

OCCURS 


COMMON CAUSES 
AND FEATURES 


1 . 


Moderate 

South 

westerly 


SE over 
mixed* 
layer 


190- 

310 


Low 


spring Seabreeze/easterly 
& fall diurnal balance, 

mesoscale eddy 


2 . 


Weak 

Westerly 


Weak 


240- 

290 


Low 


winter , 
or with 
summer 
stratus 


Seabreeze with 
mild land-sea 
temperature 
difference 


3 . 


Summer 

North 

westerly 


Weak 


290- 

340 


Low/Mid 


After 

noon 


Pacific High & 
Monsoon 


4 . 


Evening 

Northerly 


Weak 


340- 

070 


Mid 


Non- 

winter 


Veering behind 
Seabreeze front 


5. 


Stagnation 


Weak 


000- 

360 


None 

( surface) 


Near 

dawn 


Weak Seabreeze/ 
slope drainage 
balance .Usually 
not stationary 


6. 


Nocturnal 

Easterly 


varies 


040- 

160 


None 

( Surface) 


winter 


Slope drainage, 
winter monsoon, 
& diurnal 


7 . 


Synoptic 

Easterly 


Strong 

E 


040- 

160 


High or 
None 


spring 
& fall 


Warm, dry Santa 
Ana 


8 . 


Storm 

Southerly 


SW w/ 
upper 
trough 


090- 

200 


High or 
None 


winter 


Fronta 1 
passage 


9. 


Winter 

North 

westerly 


Strong 

NW 


290- 

030 


Mid/High 
or None 


After 
noon or 
storms 


Postfrontal or 
Nevada Low if 
strong, Pac High/ 
Seabreeze if weak 


10. 


Baja 


■? 


090- 


Mid 


fall 


Surface low over 



Southerly 170 Baja California 



Low = 100m 
Mid = 400m 
High = 600m 
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with the above considerations in mind we selected the following 
cases from the 1966 data period: 



Case Date 



28 


February 


16 


31 


February 


24 


48 


March 


29 


55 


April 


21 


87 


June 


13 


90 


June 


21 


91 


June 


22 


110 


June 


22 



All the flow types in Table 2.2.1, save for the Baja Southerly, can 
be viewed within the context of seasonal and diurnal variations in 
the Vandenberg area according to Table 2.2.2, also taken from the 
Handbook. Perhaps as many as five or six of the above flow types 
are represented among the eight cases. As is consistent with 
Vandenberg climatology, the vast majority of the Mt. Iron flows 
were either seabreeze or synoptic, post-frontal northwesterlies. 
These are represented in cases: 31, 55, and 90. Weaker seabreeze 
westerlies are also present in cases, 28 and 110. The evening 
northerly and perhaps partial stagnation is indicated in cases, 48 
and 91. The remaining case, 87, is a seabreeze southwesterly, 
perhaps initiated by the positioning of the edge of the local 
stratus deck rather than competition between the seabreeze and a 
synoptic easterly. 

Storm southerlies, nocturnal easterly drainage, synoptic 
easterlies, and Baja Southerlies are not indicated in the available 
data. The first of the three cases which might be construed as 
storm southerlies does not appear in the Mt. Iron data set until 
MI-175 from SLC-6 in January 1967 during Phase II. Since data 
collection during storm periods is difficult or impossible, perhaps 
this was intentionally avoided until a high level of field 
competence was achieved. The scarcity of easterly drainage flow 
cases during both Phase I and Phase II is not surprising, since 
most data were taken during daylight hours. 

The upper air data was also valuable for visualizing the temporal 
and spatial evolution of the flow over Vandenberg (e.g, sea-breeze 
front location, marine layer depth, drainage, fog or stratus 
effects) and aided in determining the flow type analysis. 

Since a number of rarer flow types characterized in the Handbook 
were not available within the data set, we also tried to skew the 
sampling to span all three available seasons and some night time 
cases. Thus, cases 31, 48, 90, and 91 were nocturnal releases. 
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while cases 28, 55, 87, and 110 were in daytime. Cases 28 and 31 
occurred during late winter, 48 and 55 in spring, and 87, 90, 91, 
and 110 occurred in early summer. Tracer was released from SLC-6 
exclusively during the fall season. 

Seven of the eight cases were releases from SLC-4 , with only case 
55 being a release from SLC-5. Though not from SLC-4, case 55 was 
chosen because it was the only case indicating strong convection. 
In Phase I there were some instances of plumes from SLC-5 largely 
confined to flow up Honda canyon during the usual mid-morning 
Seabreeze westerlies. But this pattern was not seen from SLC-4 and 
we did not include them. 



2.4 Case Flow Analysis 

The eight Mountain Iron cases extend from February to June, with 
marine boundary layer heights varying between 160 and 1300 meters. 
The following are capsule descriptions of the meteorology. 



MI-28 



Date 

Release time 
site 

duration 
Flow Type 
Inversion 



16 February 1966 
1200 local 
VIP #1 
15 min 

Seabreeze weak westerly (SWW) 
Deep, weak, at 1300 meters? 



Flow analysis : The small horizontal variation in the 

temperature profile suggests modestly higher off-coast pressures 
indicative of a typical noontime, weak winter Seabreeze. Wind 
vectors at WT301 and the Boathouse suggest that the flow diverges 
around the high terrain across Honda Ridge, where hill induced 
pressure gradients create the surface layer acceleration, evident 
at the 500m level at WT055. The veering indicated at Telemetry 
Peak and WT014 are probably due to local divergence initiated by 
the high terrain. Even at these modest wind speeds the evident 
flow distortion generated by Honda Ridge suggests a lower, stronger 
inversion than is indicated by sounding records. 



MI-31 



Date, release time 
Release time 
site 

duration 
Flow Type 
Inversion 



24 February 1966 
1845 local 
VIP #1 
15 min 

Synoptic Northwesterly (NW) 
Unknown, probably deep 
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Flow analysis : Wiresonde temperatures were available only to 

220m. Hence, the lack of upper wind or cloud data precludes 
certainty. However, the strong northwesterly flow, wind speeds and 
time of year suggest post-frontal forcing with modest veering 
within Honda Canyon and surface layer acceleration across Honda 
Ridge. The lack of flow divergence to the westward, lower portion 
of Honda Ridge also suggests a high inversion, consistent with 
post-frontal conditions. The abscence of backing once the flow has 
passed beyond the ridge is consistent with this interpretation. 



MI-48 



Date 

Release time 
site 

duration 
Flow Type 



29 March 1966 
2315 local 
VIP #1 
5 min 

Nocturnal Drainage (ND) 



Inversion : Elevated subsidence type at 400m near release 
site but with large variability. For example the inversion above 
the boathouse appears low, perhaps due to radiation fog. 



Flow analysis : Tower winds show veered, northwesterly 

remnants of a seabreeze disturbed by more veering across Honda 
Ridge. The inversion lid forces the flow \/est over the lower 
portion of Honda Ridge. South of the ridge, continuity forces the 
winds east again, as seen at the Boathouse. Downslope drainage 
flow into the canyons is seen at low levels, with complicated 
reverse flow due to radiation fog which forms in the larger 
drainage pool areas such as upper Honda Canyon. 



MI-55 



Date 

Release time 
site 

duration 
Flow Type 



21 April 1966 
1109 local 
Bldg 529 (SCOUT) 
30 min 

Northwesterly 



Inversion : Weak subsidence type at 650m. 

Flow analysis : Well-mixed boundary layer up to the inversion 

with winds maximizing at 200-300m, then decreasing aloft. Humidity 
profiles indicate cloud cover over Honda Canyon and Scout D. In 
view of the unstable lapse rates indicated by both the rawin and 
wiresondes, this suggests spring cumulus convection, consistent 
with the short plume footprint. Flow distortion across Honda 
Ridge, similar to MI-48, is still evident but weaker, perhaps due 
to the higher, weaker inversion. 
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MI-87 



Date 

Release time 
site 

duration 
Flow Type 



13 June 1966 
1310 local 
VIP #1 
30 min 

Southwesterly seabreeze 



Inversion 



Low subsidence type at 250m. 



Flow analysis : The flow is weakened and backed from its usual 
summertime pattern near the surface. However, rawinsondes and 
towers along the ridge lines show that the higher level winds 
retained the usual northwesterly seabreeze character for this time 
of day. This is the opposite pattern from that expected for 
southwesterlies induced by inland high pressures, the kind which 
give rise to hot e terly Santa Anas further south. In this case 
the seabreeze has ought a fog front with it, leading to a low 
level, thermally based, southwesterly "fog breeze". 



MI-90 



Date 

Release time 
site 

duration 
Flow Type 



21 June 1966 
2300 local 
VIP #1 
3 0 min 

Northwesterly summer seabreeze 



Inversion 



High subsidence type at 600m. 



Flow analysis : The high inversion, clear skies and strong 

winds this late into the evening suggest some frontal augmentation 
of the usual seabreeze. The usual flow distortion over and around 
Honda Ridge is again shown by the wind vectors at WT301, 055, and 
057. Tripling of the speed at WT055 suggests strong flow funneling 
near the inversion, as well as acceleration over hills, due to hill 
induced pressure gradients (Jackson and Hunt, 1979) . The anomalous 
southerly flow at the mouth of Honda Canyon is hard to explain, 
except as instrument error. However, we speculate that the high 
winds at the Boathouse may be due to a lowered inversion height and 
consequent high speed funneling south of Honda Ridge. The boundary 
layer may be lower south of the ridge because the ridge obstructs 
and dams the northwesterly flow, inducing a Scorer type of 
hydraulic jump across the ridge. Rippling of the inversion base by 
leeside gravity waves may augment this effect. 
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MI-91 



Date 

Release time 
site 

duration 
Flow Type 



22 June 1966 
0203 local 
VIP #1 
30 min 

Evening Northerly 



Inversion 



Mid-level subsidence type at 400m. 



Flow analysis : This release occurred three hours after case 
90. In the interim winds have veered from northwesterly to 
northerly and weakened somewhat, due to the lack of thermal 
forcing. This follows the summer seabreeze pattern outlined in the 
Handbook. According to rawinsondes, the veering occurred around 
0000 local time and appeared more evident above 600m, while the 
inversion base fell sharply (by "250m) . These changes may be 
consistent with post-frontal relaxation of the hydraulic jump 
mentioned in the case 90 analysis. Even so, northerly forcing 
remained intense enough to avoid any drainage or higher level 
easterly return flow tendencies, as often appear later at night. 
Again, the usual flow distortion over and around Honda Ridge shows 
up in the wind vectors at WT301, 055, and 057. And the wind vane 
at the mouth of Honda Canyon still seems in error. 



MI-110 



Date 

Release time 
site 

duration 
Flow Type 



22 June 1966 
1055 local 
VIP #1 
30 min 

Seabreeze Weak Westerly 



Inversion 



strong subsidence type at 160m 



Flow analysis : This release occurred eight hours after case 

91. In the interim winds have veered completely around to 
westerly. Inland heating has burnt back the overnight stratus 
cover to near shore, above which remains the ever present 
subsidence inversion. However, the seabreeze is still weak, due to 
the limited heating period, nor has Coriolis induced veering 
really begun yet. This follows the summer seabreeze pattern 
outlined in the Handbook. According to rawinsondes, winds backed 
rapidly above the inversion to a light southeasterly, also typical 
of the general seabreeze circulation pattern at Vandenberg. 
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LINCOM-RIMPUFF MODEL DESCRIPTION 



3 . 



3.1 LINCOM 

There are two general classes of models, diagnostic (steady state), 
and prognostic (time marching) . Diagnostic codes are suited for 
nowcasts where immediate toxic hazard corridors (THCs) are needed 
for accidental spills. Prognostic codes are more suited for 
forecasts of hypothetical spills or planned launches. This 
distinction blurs, if a diagnostic model is fast enough to provide 
frequent updates based on perhaps a five minute or better basis, or 
if faster-than-real-t ime prognostic codes are run continuously and 
frequently nudged by new input data. LINCOM is a five level, 
diagnostic boundary layer flow model, which when combined with a 
puff dispersion model such as RIMPUFF, is fast enough to allow such 
updates on a two minute basis. LINCOM 1.0 itself requires about 
one minute for a typical case at 500m resolution over a 40 x 60 km 
domain on a 386/33 PC. 

LINCOM was developed by Ib Troen at the RISO National Laboratory in 
Roskilde, Denmark (Troen, 1986). Other recent contributers have 
been George Lai, Ray Kamada , and Torben Mikkelsen. 

LINCOM uses linearized versions of the u, v, and w components of 
the momentum equation, the boundary layer mixing form of the mass 
continuity equation and the equation of state. LINCOM versions 1.0 
and 1.1 both neglect the temperature equation (also termed the 
thermodynamic energy equation) . Thus, neither stable nor unstable 
conditions are treated within the governing equations. However, 
through an objective analysis scheme applied after the solutions 
from the governing equations are obtained, LINCOM is designed to 
match the tower data exactly, with dynamic interpolation between 
towers. That is, at grid points between the tower inputs, LINCOM 
adjusts the winds by combining linearized dynamics with a terrain 
modified, inverse-distance-square based objective analysis. So 
LINCOM conforms to the actual winds and the thermodynamic and 
dynamic forcings, to the level of accuracy and resolution provided 
by the towers. Since LINCOM operates at five levels within the 
boundary layer (default values are: 6', 54', 0.2 Zi, 0.5 Zi, and 

0.8 Zi) , levels above the towers use similarity based 
extrapolations. 

LINCOM 1.0 was used for the Mt. Iron study. However, in version 
1.1, the resultant wind field is not completely anchored to the 
tower data, as is discussed below. 

LINCOM is speed optimized by two techniques. First, we transform 
the Vandenberg terrain from grid point to Fourier spectral space 
wherein the governing equations are solved. In this way the same 
resolution can be achieved using far fewer sine/cosine waves as 
grid points. Second, instead of numerically solving the governing 
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equations at each grid point, linearization lets us pre-calculate 
symbolic solutions so that only the coefficients for each wave 
number need further adjustment. Together these two techniques allow 
more than an order of magnitude speed increase over standard finite 
difference methods. 

In addition, for linearized neutral flow, Lai (unpublished) showed 
recently that the total flow can always be expressed as a linear 
combination of independent orthogonal components of the mean and 
perturbed flows along just the x and y axes. Since the orthogonal 
flows are pre-calculable , we can simply find the mean wind speed 
and direction which best matches the Vandenberg tower winds and 
apply them to pre-calculated solutions. In version 1.0, a 
"look-up" mode employs a set of 72 wind fields spaced in three 
degree increments around 360® of arc which are interpolated to 
obtain the final field. However, in version 1.1, the best linear 
combination of two pre-calculated orthogonal fields is provided. 

Trade-offs are involved in formulating such a model. Since LINCOM 
does not account for time tendencies in the velocity, 
dynamical forcing due to terrain is effectively propagated 
instantaneously over the entire domain, rather than at some 
realistic set of phase speeds. Linearization of the governing 
equations also means that second and higher order perturbation 
terms are neglected. In a prognostic scheme, such solutions 
would gradually diverge from reality. However, for a diagnostic 
model, where in effect only the present time step is being 
considered, this issue is not very significant. 

Another aspect of linearization in LINCOM is that for each grid 
point, the vertical velocity is basically assumed equal to the 
horizontal velocity at that site times the sine of the terrain 
slope. This approximation only holds well for slopes of less than 
twenty degrees. Thus, the aspect ratio of the terrain (typical 
height/typical length) must be relatively small compared to unity. 
ACTA's analysis of the Vandenberg terrain concludes that at 500 
meter resolution, hardly any of the terrain shows slopes greater 
than twenty degrees (Conley et al., 1990). Thus, the quantitative 
limit for the vertical velocity assumption is seldom exceeded at 
Vandenberg . 

Transformation of the domain to spectral space involves the 
assumption of periodic boundary conditions. So wave energy 
leaving one boundary will immediately re-enter the domain from the 
opposite boundary. For this reason, a 10 km buffer zone was 
created around the Vandenberg domain in which the terrain is 
gradually relaxed on all sides to sea level. This has the effect 
of greatly damping such wave motions. However, the effect of 
periodic boundary conditions cannot be avoided completely and we do 
not have a quantitative assessment of its influence on the flow as 
yet . 
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Again, the effect of the objective analysis in LINCOM 1.0 is to 
anchor the resultant wind field to the wind vectors given by the 
Vandenberg base tower system. In LINCOM 1.0 Monin-Obukhov 
similarity theory is used to extend the tower vectors upward into 
the outer boundary layer. In effect this turns LINCOM into a 
dynamically based interpolation/extrapolation scheme based on the 
tower winds. However, the tower winds may be influenced by effects 
which are more local than can be resolved by the implemented 500 
meter grid spacing, such as buildings and local terrain. This 
constitutes a type of aliasing. Since the objective analysis is 
not inherently mass conserving, such local streamline divergence 
indicated by the towers can cause departures from the continuity 
constraint implemented in the model physics. However, since LINCOM 
1.1 is a best fit to, rather than an anchoring by the tower 
vectors, strict mass conservation is retained, while avoiding the 
aliasing issue. 

Most of the available towers are clustered near the coast on 
Vandenberg property. Other towers are available from the Santa 
Barbara Air Quality Management District's set of towers, some of 
which are actually maintained by the regional oil refineries. 
There are also two local buoys maintained by NOAA and towers atop 
a few of the local oil drilling platforms. However, these sources 
have not been incorporated as yet into the on-line data stream. 



3 . 2 RIMPUFF 

There are four classes of plume diffusion models which are being 
considered for use at Vandenberg. In ascending order of 
computational demands they are: 1) steady state plume model with 
gaussian lateral concentration profile, 2) mixed layer scaling 
model, 3) lagrangian puff model with gaussian radial concentration 
profile, 4) lagrangian particle model with particle motions 
governed by a mean wind plus a random turbulent component with a 
gaussian velocity profile. RIMPUFF is a lagrangian puff model. 

RIMPUFF is structured to handle multiple simultaneous sources. 
Release points can be located anywhere within the 3-D grid and can 
be specified individual release rates, release times and heat 
production. A series of puffs is released to simulate a 
continuous plume. At each time step, a book-keeping algorithm 
tracks the advection, diffusion, and fractional deposition of 
each puff in accordance with local meteorological parameters. 

The model calculates the concentration at each grid point by 
summing the contributions from all the surrounding puffs. These 
concentrations can simply be updated or also accumulated to provide 
dosage. The model output consists of individual puff locations and 
grid concentrations at user specified time intervals. From these 
results, a graphics program produces puff plots and concentration 
or dose isopleths. 
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RIMPUFF employs the LINCOM mean flow field to advect contaminant 
puffs downstream. Within-puff growth is controlled by local 
turbulence level, using two-particle relative diffusion theory. 
Local turbulence intensity depends on tower based climatological/ 
directional parameterizations of Oq, the standard deviation of wind 
direction, for lateral turbulence, and Pasquill-Gif ford stability 
classes modified for complex terrain for the vertical turbulence. 

RIMPUFF treats plume bifurcation in complex terrain by letting 
puffs split when they exceed the size of the grid spacing. I.e., 
for Vandenberg an initial 100m puff will grow to the 500m LINCOM 
grid spacing before it splits horizontally into five 250m puffs 
(pentaf urcation) . Vertical trifurcation can occur independently. 
The practical splitting limit on PC based computers is several 
hundred puff progeny. Mass, momentum, and center concentration (of 
the mother puff) are conserved in the splitting process. If the 
mean flow is not frequently updated and remains static, RIMPUFF 
adds stochastic (Langevin based) advection to simulate puff meander 
due to sub-grid turbulence. This is particularly important when 
puff siblings are closely spaced shortly after splitting. 

Since the series of puffs is meant to simulate a continuous plume, 
the effective stack height after plume rise is taken from standard 
plume rise models for hot spill cases. 100% reflection is assumed 
at a user specified inversion height which parallels the terrain. 
The final rise height for each puff is a function of the 
atmospheric stability and windspeed at the time and height of the 
release as given by LINCOM. 

For Mt. Iron the ground level reflection coefficient was set to 
100%. Dry deposition is calculated using the source depletion 
concept according to stability and windspeed. Wet deposition 
accounts for rain intensity and duration in time and space, using 
a 1/r^ objective analysis, where r is the distance from the grid 
point to the measurement station. 

As a puff model RIMPUFF has advantages over gaussian plume and 
similarity scaling models. It can handle non-homogeneous terrain 
during non-stat ionary conditions. As plumes extend to farther 
distances and longer travel times, they are more likely to 
encounter and respond to such variations. In complex terrain, 
non-gaussian profiles can occur which RIMPUFF can treat, such as 
bifurcated plumes, pooling, non-homogeneous channeling, and 
fanning . 

Gaussian plume and mixed layer scaling models must also assume some 
kind of profile for the vertical concentration distribution. A 
gaussian vertical profile may be modified to account for surface 
reflection. This is appropriate for short distances before a 
significant amount of reflection occurs from an upper level 
inversion. At long distances a uniform vertical concentration may 
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be assumed. In fact, neither of these simple assumptions is 
entirely reasonable, since plumes will spread vertically until an 
inversion is reached, then reflect and mix downward in a gradual 
and complex fashion. Non-stationary lagrangian puffs and particles 
should be able to model such complex phenomena more accurately. 

Puff model accuracy is mitigated, especially near the source, by 
finite puff size and assumed radially gaussian concentration 
profiles. Lagrangian particle models improve on these aspects. But 
they reguire a much larger particle swarm and more computer time. 
RIMPUFF's stochastic puff advection scheme is also not well tested 
yet. On the other hand, lagrangian particle models assume a 
gaussian velocity component to diffuse the particles. Yet, the 
vertical velocity distribution has positive skewness and therefore 
is not gaussian under convective conditions. Baerentsen and 
Berkowicz (1984) and Misra (1982), however, claim better results 
with an adjustable double gaussian distribution for their 
lagrangian models. 



4.0 DATA/MODEL COMPARISON METHODS 



4.1 Initialization 

With a model to model comparison, it is trivial to create output 
grids with compatible domain size and grid spacing to 
guantitatively compare the results. However, a model to data 
comparison such as Mt. Iron requires some manipulation before 
quantitative results may be obtained. 

Among the products of the Mt. Iron Experiment was a set of hand- 
drawn isopleths defining measured dosages for each case, given in 
terms of time-integrated concentrations normalized by release rate 
(units of lO^sec/m^). However, these results were not immediately 
compatible with the format required for comparison with LINCOM/ 
RIMPUFF which defines dosages at distinct points on a grid. Thus, 
the eight cases from Mt. Iron were digitized onto a 119 x 84 evenly 
spaced grid at 100m mesh spacing, and formatted for the commercial 
SURFER PC plotting package (Golden Software, 1989) . This domain 
size was large enough to encompass the Mt. Iron and LINCOM/RIMPUFF 
plumes yet small enough to comply with the grid size constraints of 
SURFER (see figs. 1 and 2) . 

Inspection revealed that some of the hand-drawn dose isopleths from 
Mt. Iron were extended beyond the scope of the bag sampling 
network, presumably to ensure closure of the contour lines. Such 
extrapolation beyond the measurements in cases: 28, 48, 87, 90, and 
110 cannot really be justified, since it relies on subjective 
judgement. In these cases both the modeled and measured isopleths 
were truncated along the last measurement line. For instance, in 
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case 28, the measured plume passed over the sampling line along 
Arguello Road, but did not reach Santa Ynez Road. Since the exact 
loci for dose isopleths between these two sampling lines was not 
known, the isopleths were truncated at Arguello Road for the DICA 
and other comparisons. 

For the following cases the cut-off lines were: 



Iron Case 


Plume Cutoff 


28 


Arguello 


48 


Honda Ridge 


87 


Arguello / Santa Ynez 


90 


Tranquillon / Honda Ridge 


110 


Arguello 



Measured dosages along the Coast Road in the vicinity of Sudden 
Ranch for case 91 were sufficient to resolve the plume, whereas for 
case 48, measured dosages were too small and the plume was 
truncated at Honda Ridge. 



4.2 DICA and Other Merit Scoring Methods 

The Dose Isopleth Correlation Area (DICA) scoring method is a 
conceptually simple performance measure we created to compare the 
coincidence between LINCOM/RIMPUFF dose isopleths and the hand 
drawn ones, based on Mt. Iron. From simple set theory, for each 
case and each dose isopleth within a case, DICA divides the total 
area of intersection between each modeled and measured isopleth by 
the total area of union. However, before being summed, it weights 
each intersect by the measured to modeled isopleth value ratio. 
Thus, if there are k number of isopleths with modeled and measured 
doses, aj, aj , and modeled and measured areas, Aj, Aj, the DICA 
algorithm is 



k k k k 

Z E (a /a )A n A + Z Z (a /a )A fl A 

j=l i=l i j i j i=l j=l j i j i 

a <a i=j 

i j a <a 

j i 



(4.2.1) 



k k 

Z A U Z A 
i=l i j=l j 
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This yields a single value, ranging between zero and one, for the 
amount of overlap between modeled and measured results. Perfect 
congruence yields unity, while no overlap yields zero. 

DICA presents a fairly stiff test of "goodness of fit". That is, 
the congruence must be guite good to obtain DICA values above 0.5. 
Thus, we suggest the following interpretation of DICA scores, 



DICA 


VALUE 


FITTING ACCURACY 


0.5 


- 1 


excellent 


0.2 


- 0.5 


good 


0.05 


- 0.2 


fair 


0.00 


- 0.05 


poor . 



In addition to DICA two other performance measures were utilized as 
suggested by Hanna (1991), fractional bias. 



FB = (Do - D^)/(Dq + D,J , (4.2.2) 

and normalized mean square error. 



NMSE = (Do - D^)/DoD,„ , (4.2.3) 



where Do and D^ are the modeled and measured dosages at each grid 
point. Programs were written to determine DICA, FB, and NMSE 
values for each of the eight cases. 

FB is useful for showing which plume predicts the largest dosages. 
A zero FB indicates that the model predicts on average the same 
dosage as was measured. The NMSE should be behaviorally similar to 
DICA, except that lower NMSE should imply better model predictions, 
while higher DICA values indicate a better fit (see discussion of 
relative merits in section 5.2). 

All values which exceeded the lowest isopleth value of 20 are 
included. If either the RIMPUFF or Mt. Iron dosage value at a 
given grid point exceeded the minimum threshold value, then that 
grid point was included in the averaging process for the FB and the 
NMSE scores. 

We also show a standard comparison of observed versus modeled 
centerline dose exposure levels, including factor of two and factor 
of four error limits (fig. 27). 
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5. DATA/MODEL COMPARISON RESULTS 



5.1 LINCOM Winds, RIMPUFF/Data Isopleths, DICA and Other Results 



Results of each Mt. Iron case and LINCOM/RIMPUFF simulation are 
plotted in figs. 3-27. Dose isopleths are given in terms of time- 
integrated concentrations normalized by release rate (10^ sec/m^) . 
Locations are given in Universal Transverse Meridian System (UTMS) 
coordinates at 100 m resolution. The release location for case 55 
was at the mouth of Honda Canyon near SLC-5. The other seven 
releases were from VIP #1 at the northwest corner of the SLC-4 area 
along the Coast Road. 

The DICA, FB, and NMSE values are presented in Table 5.1 for the 
eight cases. 



Table 5.1.1 



Comparison of 



model 



performance measures. 



Case 


Number 

of 

points 


Do 


1 

I 


(Do-Dp)- 


FB 


NMSE 


DICA 


28 


575 


420 


684 


1, 175,025 


-0. 477 


4 . 082 


0. 369 


31 


1038 


231 


181? 


125,231 


+0.207 


2.863 


0.266 


48 


1499 


171 


225 


135,991 


-0.233 


3.391 


0.268 


55 


692 


132 


142 


133 , 124 


-0 . 069 


7 . 063 


0. 121 


87 


692 


264 


355 


265,698 


-0.294 


2 . 825 


0.535 


90 


1589 


285 


460 


837,967 


-0.470 


6.375 


0.328 


91 


2209 


177 


187 


387 , 625 


-0.054 


11.700 


0.369 


110 


517 


303 


470 


484,409 


-0.432 


3 . 395 


0.472 1 
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Fig. 3 Mt . Iron Cbs© 28 plum© isopl©ths sndl tow©ir wind v©ctoirs 
from th© Mt. Iron R©port Phas© I Vol. I. 
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Fig. 4 LINCOM Case 
shown at every fourth 



28 (weak seabreeze westerly) wind vectors 
grid point (2km). Domain is 25 x 40 km. 
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LINCOM vs. Mi. Iron Case 28 
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Fig. 5 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 28 at 
100m resolution. Domain is 8.4 x 11.9 km. 
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Fig. 6 Mt. Iron Case 31 plume isopleths and tower wind 
from the Mt. Iron Report Phase I Vol. I. 
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Fig. 7 LINCOM Case 31 (synoptic northwesterly) wind vectors shown 
at every fourth grid point (2km) . Domain is 25 x 40 km. 
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LINCOM vs. Mi. Iron Case 31 
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Krtx, UTMS 

Fig. 8 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 31 at 
100m resolution. Domain is 8.4 x 11.9 km. 
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Fig. 9 Mt . Iron Case 48 plume isopleths and tower wind vectors 
from the Mt. Iron Report Phase I Vol. I. 
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Fig. 10 LINCOM Case 48 (nocturnal drainage) wind vectors shown at 
every fourth grid point (2km) . Domain is 25 x 40 km. 
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LINCOM vs. Mt. Iron Case 48 
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Fig. 11 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 48 
100m resolution. Domain is 8.4 x 11.9 km. 
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Fig. 12 Mt. Iron Case 55 plume isopleths and tower wind vectors 
from the Mt. Iron Report Phase I Vol. I. 
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LINCOM Case 55 (convective seabreeze northwesterly) wind 
shown at every fourth grid point (2km) . Domain is 25 x 40 
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LIN COM vs. ML Iron Case 55 

717 718 719 720 721 722 723 724 725 




3838 

3837 

3836 

3835 

3834 

3833 

3832 

3831 

3830 

3829 

3828 

3827 



14 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 55 at 
resolution. Domain is 8.4 x 11.9 km. 




Fig. 15 Mt. Iron Case 87 plume isopleths and tower wind vectors 
from the Mt. Iron Report Phase I Vol. I. 
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Fig. 16 LINCOM Case 
shown at every fourth 



87 (fog breeze southwesterly) wind vectors 
grid point (2km). Domain is 25 x 40 km. 
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LINCOM vs. Ml. Iron Case 87 
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Fig. 17 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 87 at 
lOOtn resolution. Domain is 8.4 x 11.9 km. 
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Fig. 18 Mt. Iron Case 90 plume isopleths and tower wind vectors 
from the Mt. Iron Report Phase I Vol. I. 
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Fig. 19 LINCOM Case 90 (seabreeze northwesterly) wind vectors 
shown at every fourth grid point (2km). Domain is 25 x 40 km. 
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LINCOM vs. Mt. Iron Case 90 
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Fig. 20 LINCOM/ RIMPUFF and Mt . Iron dose isopleths for Case 90 
100m resolution. Domain is 8.4 x 11.9 km. 
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Iron Report Phase I Vol. I. 



vectors 



40 



0 5 10 15 20 25 




Scale 

— ^ 11.00 m/s 



Fig. 22 LINCOM Case 91 (evening northerly) wind vectors shown at 
every fourth grid point (2km) . Domain is 25 x 40 km. 
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LINCOM vs. Mi. Iron Case 91 
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Fig. 23 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 91 at 
100m resolution. Domain is 8.4 x 11.9 km. 
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Fig. 24 Mt. Iron Case 110 plume isopleths and tower wind vectors 
from the Mt. Iron Report Phase I Vol. I. 
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Fig. 25 LINCOM Case 110 (seabreeze weak westerly) wind vectors 
shown at every fourth grid point (2kin). Domain is 25 x 40 km. 
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LINCOM vs. Mt. Iron Case 110 
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Fig. 26 LINCOM/RIMPUFF and Mt. Iron dose isopleths for Case 110 at 
100m resolution. Domain is 8.4 x 11.9 km. 
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Mi. Iron vs. RIMPUFF Dosages 




Observed Mt Iron Doso.ge (s/m**3) 



Fig. 27 Scatter plot of normalized centerline dose exposures for 
Mt. Iron versus LINCOM/RIMPUFF . Units are sec/m- 
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5.2 Discussion of Results 



As shown in figs. 3 - 18, the LINCOM generated wind fields conform 
quite closely to those observed. This level of conformity is not 
usually seen in diagnostic models and is mainly due to the "tower 
anchoring" feature discussed in section 3.1. Since the towers were 
not located at exact grid points and tower influence decays with a 
modified inverse distance squared, the more exaggerated flow 
divergence indicated by the towers is mitigated on the gridded wind 
field by LINCOM dynamics. However, in cases 90 and 91 this 
conformity becomes overzealous, since Tower 012 at the mouth of 
Honda Canyon was probably in error (figs. 13, 15). 

As an example of the "aliasing" problem discussed in section 3.1, 
LINCOM is also overly influenced by towers at the bottoms of 
canyons where the winds are very locally controlled but do not 
conform to the bulk of the boundary layer flow. The imposition of 
strict mass conservation as in LINCOM 1.1 should retain the main 
flow features but would not resolve the local flows. In fact, 
tower observations will generally suggest more flow divergence 
than is displayed in modeled wind fields because the tower winds 
respond to all forcings, including the ones that are of smaller 
scale than the model grid resolution such as turbulence and very 
local terrain. 

This leads to a partial dilemma for tracer studies in hilly, 
complex terrain. As at Vandenberg, available roads for fixed or 
mobile samplers are often along the bottoms of canyons. If the 
canyon transects the plume, canyon samples often will show wider 
plume widths, reflecting the tendency for local along-canyon flow 
(see LVDE, 1991) . However, most of the plume will be relatively 
unaffected and float over such canyons. Without very high 
resolution and complete prognostic dynamics, numerical models may 
not agree well with data when the sampling does not conform to the 
bulk of the flow. Of course aircraft data would tend to mitigate 
this sampling flaw. 

The height of the inversion base is another important influence on 
the flow field. LINCOM/RIMPUFF assumes that the inversion base 
parallels the terrain. For daytime unstable conditions, this is a 
reasonable assumption as shown in Kamada et al (1990). However, 
local inversion height changes are suggested in cases 48 and 90, 
(see section 2.4). A lowered inversion may be responsible for the 
high winds seen at the Boathouse in case 90 (fig. 13). The 
inversion base may also be somewhat flatter than terrain parallel 
under neutral or mildly stable conditions. These features require 
a complete physics package and cannot be resolved by the dynamics 
of simpler diagnostic models. On the other hand, in this case the 
"tower anchoring" feature in LINCOM assists modeling accuracy in 
the winds south of Honda Ridge, in spite of the crude inversion 
height assumption. 
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Figs. 19 - 26 and the corresponding DICA, FB, and NMSE scores 
suggest that LINCOM/RIMPUFF matches the Mt. Iron isopleths fairly 
well, except for case 55 where the DICA score falls below 0.2. 
However, the fractional biases indicate that the LINCOM/RIMPUFF 
plumes are longer than Mt. Iron in seven of the eight cases. I.e., 
doses are higher than those measured at longer downwind ranges, 
particularly cases 28, 90, and 110, where FB are less than -0.4, 
even though downwind ranges for the higher dose isopleths are quite 
comparable. This is in keeping with the AFTOX/Mt. Iron and 
WADOCT/Mt. Iron comparisons of plume centerline dosages (Kunkel and 
Izumi, 1990; Kunkel, 1991). It is also consistent with the LVDE 
results . 

As in the above reports, we suggest that the zinc sulfide wet 
aerosol tracer was not physically inert as previously assumed, but 
probably was gradually adsorbed into the ground canopy as the plume 
was advected downwind. This is consistent with oral reports that by 
the conclusion of the Mt. Iron releases, pools of dried tracer 
material surrounded the release sites. 

For most cases, the DICA and NMSE scores yield parallel results. 
For instance, case 87 shows the best fit with both methods. 
However, the DICA and NMSE scores diverge for cases 31, 90, and 
particularly 91. The NMSE was much higher for case 91 than the 
other cases, even though the DICA and FB scores were fairly good. 
As in 31 and 90, the NMSE in case 91 involves a large number of 
points (column 2 of Table 5.1.1). However, both the measured and 
simulated average dosages are low (columns 3 and 4), as is the 
total error (column 5) . But the NMSE is large here because the dose 
average is so small. Case 90, with a much larger average dosage 
and error, has a much smaller NMSE. This suggests that normalizing 
the error by the average isopleth value may skew the results to 
favor cases with large average dosages, despite having large errors 
prior to normalization. 

On the other hand, the DICA algorithm does not seem to suffer from 
this problem. It sorts the LINCOM/RIMPUFF model results for the 
eight cases in much the same order as derived by a qualitative 
visual inspection, but more accurately. For example, a visual 
inspection (Fig. 23) suggests that case 87 shows the best fit. It 
also has the highest DICA score; case 55 (fig. 22) shows the worst 
fit and lowest DICA score. However, it is difficult to distinguish 
case 87 from 110 or to compare two rather different cases, such as 
28 and 91, for congruence by eye, while DICA does so easily. Given 
these strengths, we suggest that DICA be utilized as a standard 
performance measure where applicable. 

Save for cases 55 and 90 (figs. 22 and 24), the modeled plume 
thicknesses are comparable to the Mt. Iron estimates. In cases 55 
and 90 RIMPUFF overestimates the lateral dose footprint. 55 is a 
special case discussed in more detail below. However, we suspect 
that the problem in case 90 stems from underestimating the vertical 
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diffusion. Wind/temperature vertical profiles for case 90 suggest 
near neutral conditions with significant mechanical turbulence. 
For such conditions, RIMPUFF's vertical turbulence parameterization 
will tend to underestimate the extent of the shear based 
turbulence. In general, it is also hard to assure accuracy for 
tower based temperature profiles on a routine operational basis, 
due to the need for frequent, precise sensor calibrations. Perhaps 
the surface radiation (rather than profile based) stability 
algorithm in LINCOM 1.1 should be used in RIMPUFF as well. 

Whereas the releases were nearly point sources, in each case the 
modeled plume head appears thicker, due to the initial 100 meter 
puff diameter. A smaller initial puff diameter could have been 
used. So this is not a wholly necessary model artifact. Part of 
the extra thickness is also due to the interpolation routine used 
in the plotting program, as seen from the increased thickness in 
interpolating the hand drawn Mt. Iron plots. However, the higher 
valued dose isopleths are generally broader in RIMPUFF than the 
hand drawn isopleths. This is not as true of the lower valued dose 
isopleths, since diffusion decays exponentially with distance. 

V'Je also see that whenever the plume passes over Honda Canyon, 
substantial deviations between the modeled and hand drawn isopleths 
can result. The modeled winds are influenced by towers located at 
the bottom of Honda Canyon, while in cases 31, 48, 90, and 91, much 
of the real plume seems to have floated above and beyond the 
sampling lines at the bottom of the canyon. The LINCOM/ RIMPUFF 
tandem cannot reproduce the local minimum dosage values 
seen in these cases at the bottom of Honda Canyon. Kunkel (1990) 
reports a similar modeling constraint for the WADOCT model. 
Without non-hydrostatic capability, i.e., both buoyancy and 
acceleration in the vertical momentum equation, neutral LINCOM's 
ability to simulate the vertical divergence indicated by these 
cases is quite limited. For example, the hydrostatic version of 
HOTMAC/RAPTAD also did not reproduce the local minimums seen in 
cases 90 and 91, even though HOTMAC does include buoyancy in its 
vertical momentum equation (Yamada and Bunker, 1991). It may be 
that resolutions higher than 500m are also required to show such 
effects . 

On the other hand, cases 31, 48, and 90 (figs. 20, 21, and 24) 
probably indicate aliasing in LINCOM 1.0. Presumably, the mass 
conservation integral to LINCOM 1.1 and the inclusion of buoyancy 
in LINCOM 2.0 should improve the simulation in such cases. 

Perhaps roughly 80% of the accuracy in plume dispersion estimates 
depends on accurate initial wind directions. For example, case 28 
shows that a few degrees error (or perhaps real wind shift) in 
specifying wind direction can lead to substantially lower DICA and 
higher NMSE values. 

The effect of downwind towers cannot be ascertained for the weak 
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westerly and southwesterly Seabreeze cases: 28, 87, and 110 (figs. 
19, 23, and 26), due to the lack of data beyond the sampling line 
at Arguello Rd. and consequent truncation of the plumes. 

For case 55 (fig. 22) the LINCOM/RIMPUFF plume is much longer than 
the measured one. As discussed in the case flow analysis of 
section 2.3, a convective cumulus condition with vigorous updrafts 
is indicated by the unstable lapse rates, high mid-April inversion, 
and patchy cloud cover over neighboring sites in Honda Canyon and 
Scout D. Thus, the plume may simply have lofted and dispersed over 
unsampled areas. Or a sustained local downdraft may have forced 
adsorption into the canopy. The former interpretation is probably 
more likely. 

There are data problems which render the scatter diagram in fig. 27 
less meaningful than is immediately apparent. 1) Since we compare 
only eight cases, there are only twenty-two valid data points. 2) 
Only the hand-drawn isopleths were available, rather than actual 
bag sample exposure levels. Hence, non-linear interpolation was 
needed to obtain actual values. 3) The sampling lines for Mt. Iron 
were along accessible roads rather than concentric circles. Thus, 
it was sometimes difficult or impossible to compare modeled and 
measured values when the centerlines diverged in direction. 4) 
Figure 19 seems to coincide with the fractional bias results which 
show that the RIMPUFF predictions are higher than those observed at 
lower dosage levels; actually, the predicted centerline dosages are 
higher mainly only in Honda Canyon. 



6. SUMMARY AND CONCLUSIONS 

Results from the LINCOM/RIMPUFF dispersion model were compared with 
eight representative cases from the Mt. Iron tracer study conducted 
at South Vandenberg AFB during 1966. We conclude that: 

1) The LINCOM wind fields conform fairly closely to observations 
due to "tower anchoring". Tower anchoring allows LINCOM to retain 
complex flow features which would otherwise require the use of a 
higher resolution, non-hydrostatic prognostic flow model with a 
complete physics package. However, this feature is a mixed 
blessing when local towers misrepresent the main flow. For tracer 
studies, this points to the need for aircraft sampling in addition 
to ground sampling to determine local as well as bulk flow 
features . 

2) Plume simulations based on LINCOM/RIMPUFF compare fairly well 
with Mt. Iron, but the modeled plumes are generally somewhat 
longer, as is consistent with other published reports. Perhaps 
part of this is due to some adsorption of the zinc sulfide wet 
aerosol onto the vegetation canopy during transit. 

3) Four indicators were used to compare the modeled dispersion 
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pattern with observations. The most useful appear to be scores 
from the dose isopleth correlation area (DICA) method and the 
fractional bias (FB) . We also used the normalized mean square 
error (NMSE) and a centerline scatterplot comparison. The NMSE 
seems to suffer from a tendency to overweight errors when dosages 
are relatively low. The centerline scatterplot involves too few 
data points and suffers the effects of some other Mt. Iron data 
idiosyncrasies discussed in section 5.2. 

4) Lateral plume thicknesses compare fairly well with those 
observed, except perhaps when turbulence intensities were mis- 
calculated, due to observational error, limited spatial resolution, 
or deficiencies in the turbulence parameterization algorithms. 

5) The low dosages recorded in Honda Canyon indicate that plumes 
tend to float over the canyon rather than parallel the terrain 
contours. LINCOM/RIMPUFF ' s ability to handle this feature is 
rather limited, since a complete vertical momentum equation is not 
included. Other diagnostic and prognostic model simulations of Mt. 
Iron have also shown an inability to reproduce such local dosage 
minima. Resolutions higher than 500m may be helpful in refining 
such simulations. 
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8. APPENDICES 



8.1 LINCOM Model Theory 



8.1.1 Introduction 

We simulated plume dispersion for each of the above ten flow types 
with a combined flow/puff model. LINCOM, the flow model, extends 
linear neutral potential flow theory (over a single hill) to a more 
general case involving mesoscale complex terrain. LINCOM is a 
diagnostic code which employs highly efficient solutions, which are 
spectrally based in the horizontal plane and analytic in the 
vertical. A "look-up table" of pre-computed results streamlines 
this procedure even further. Based on data from ten or more 

towers, LINCOM 1.0 can output a 100 x 140 gridded wind field for 
five vertical levels in the boundary layer at 500m horizontal 
resolution in less than one minute on a 386/33 PC. 

The flow model is initialized with mean uniform velocity, pressure, 
and potential temperature fields, (Uq, Vq, Wq, and Pq) » using data 
from each tower one at a time. Thus, at each tower the 
perturbation quantities u, v, w, and p induced by the terrain are 
all zero. LINCOM then solves the governing hydrodynamic equations 
solely in terms of the perturbation quantities; and adds the 
perturbation fields to the mean fields to obtain the final results: 
U = Uq +u, V = Vq + V, etc. So at any distance away from that tower 
the u, V, w, and p can assume non-zero values. We repeat this 
procedure for all n towers to obtain n separate initial 
perturbation wind fields. These fields are combined with a modified 
1/r* weighting scheme which warps each tower's circle of influence 
into ovals that are essentially congruent with the elevation 
contours for the surrounding terrain. This method ensures that the 
final field will exactly match the tower data, while tower 
influences extend strongly only along similar terrain elevations. 
In this way, even under differing stability conditions, a 
sufficient density of towers will tend to restore some of the 
non-neutral physics and force the final field to closely resemble 
the measured flow. 

One problem with this procedure is that mass continuity can be 
compromised. Thus, in LINCOM 1.1 we instead find the set of mean 
fields which provides the best fit to the existing set of tower 
data when added to the resultant perturbation fields. 



8.1.2 Governing equations 

To discuss the modeling ideas in more detail, we begin with a 
truncation of the full atmospheric equation set which neglects 
density variations due to effects other than gravity. This 
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standard Boussinesq set forms the basis for almost all atmospheric 
modeling. Included are the three dimensional momentum equations 
(neglecting horizontal shear) : 
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the thermodynamic energy (or potential temperature) equation, 

3e ae ae ae i ae 

— + U — + V — + W — + wr ; (8.1.2) 

at dx ay az pc az 

p 

the shallow convection form of the mass conservation or 
incompressible fluid continuity equation (applicable to boundary 
layer models) , and the ideal gas equation of state, 

au av aw 

— + — + — =0 ; p= pRT . (8.1.3) 

ax ay az 



Here, f, the Coriolis parameter is defined as 2nsin 0, where fl is 
the earth's rate of rotation. T, the potential temperature lapse 
rate, is dQ/dz. R and Cp are the gas and specific heat constants 
for air, and the Ty represent eddy stresses. We wish to simplify 
these equations because their complete solution at all relevant 
scales of motion in what is termed a direct simulation is beyond 
the forseeable capacity of computers. That is, the current direct 
simulation limit on CRAY class supercomputers is " 10^ grid points 
for mesoscale flows. A 3-D direct simulation would require " 10^^ 
grid points to resolve the larger mesoscale (lO^m) circulations, 
the dissipation of turbulent kinetic energy which occurs down to 
lO'^m, and the intervening scales of turbulent and forced motions. 
Since operational models use less than 10^ grid points, sub-grid 
scale effects must be parameterized rather than modeled directly. 
Much of the energy is resolved on the grid itself; however, 
sub-grid motions cannot be neglected because the non-linear 
eddy stress dynamically links all scales of turbulence in an energy 
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cascade. Thus, without dissipation into internal energy, that is, 
heat, the modeled turbulence kinetic energy would simply grow 
without bound. 

Other difficulties abound. For example, due to the continuum of 
turbulence scales, simple analytic forms for eddy stress do not 
exist. The better approximations are spatially dependent. They 
also involve higher order correlations between both supra and 
sub-grid velocity and pressure fluctuations whose magnitudes are 
still being debated. Non-linear advection also creates annoying 
cross terms which are usually simply ignored when the equations are 
written in non-orthogonal terrain following coordinates. Thus, 
even with supercomputers, the full or "primitive" equation 
mesoscale models are still too unwieldy for operational purposes. 

Thus, to simplify in a way which hopefully blends suitability with 
convenience, we can make linearizing approximations of the sort. 



u3u/3x = (Uq + u)d(U(j + u)/3x ” UgSu/dx 



(8.1.4) 



Fourier transforming the linearized equations to spectral space 
will then allow symbolic matrix inversion solutions for each 
wavenumber for the entire domain at once, rather than grid point by 
grid point as in finite differencing methods. However, we must 
remain aware that the linearization procedure neglects the higher 
order products of perturbed quantities, such as udu/dx. This is 
only valid if u3u/3x << Uq3u/3x. 

To suit this level of modeling we match the above approximation 
with a first order approximation of the eddy stress, i.e., = -K 

du/dz, where K is the horizontally diffusivity coefficient. 

Moreover, for stationary flows or at least if 3U/3t << U3U/3x, we 
can also neglect time derivatives. This inequality implies that 
1/T << u/L, where T and L are characteristic time and length scales 
for flow changes. For a typical T ~ 6hr and u ~ Im/s, we need an 
L << 20 km. For u ~ lOm/s, we need an L << 200 km. These are both 
reasonable for Vandenberg. 

In LINCOM 1.0 and 1 . 1 we also neglect buoyancy forces caused by an 
uneven density/ temperature field. So neutral LINCOM cannot 
diagnose thermally induced slope flows or seabreezes, except as 
indicated within the mean flow given by the towers or the later 
tower based objective analysis. this neglect allows a 

decomposition of the flow field into orthogonal components which 
are pre-calculable for a given domain and thus lead almost 
immediately to a final wind field. Also, thermal forcing at scales 
less than 10 km is often subordinate to the dynamic pressure forces 
which neutral LINCOM does account for. 
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At any rate, we are left with the following simplified, steady 
state, linearized equation set for the perturbed part of the flow: 
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continuity and 


ideal gas equations remain 


unchanged. The 1/p 



density factor has been absorbed into the perturbation pressure, p. 
As discussed, we have dropped the buoyancy term, g50/0. We have 
also dropped the small vertical terms. The solution set is 
analytic for any choice of the perturbation quantities. (u, v, w, 
p, and K) 



8.1.3 Solution method 

To solve eqns. (8.1.2, 3, and 5) we drop the Coriolis term and 
begin with a simplified 2-D version of the model in order to sketch 
the main ideas. Later, we add more complex refinements included in 
the real 3-D model. First, we Fourier transform the x component of 
velocity into wave number (k) space in the fashion. 



00 

u(k,z) = j u(x,z) exp(ikx) dx , (8.1.6a) 

-00 



1 00 

u(x,z) = J u(k,z) exp(-ikx) 3k , (8.1.6b) 

2n -00 



where eqn. (8.1.6b) is the back transform of u(k,z). Note that the 
z direction does not participate in the transform and is written 
only to remind us of the height dependence of u. Then by 
differentiation of eqn. (8.1.6a), 
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5u 



iku 



(8.1.7a) 



00 

— = ikj u(k,z) exp(-ikx) 3k = 

3x -00 



Thus, the momentum equation for u may be written as, 



00 00 

Upikj u (k , z ) exp (-ikx) 9k + ikJ p(k, z) exp(-ikx) 9k - 
— 00 -00 



9 00 

K — j u (k , z ) exp ( -ikx) 9k = 0. (8.1.7b) 

9z -00 

So the sum of the wave number components must total to zero. For 
each wave number, the 2-D version of the governing equations will 
then look like. 



ikUoU(k,z) + ikp(k,z) - K 9^u/9z^ = 0 , (8.1.8a) 

ikUQw(k,z) + 9p(k,z)/9z - K 9^w/9z^ = 0 , (8.1.8b) 



iku(k,z) + 9w(k,z)/9z = 0 . (8.1.8c) 



By virtue of the Fourier transform, the sum of the perturbation 
quantities at each wave number, k, summed over all k, is equivalent 
to their determination for each grid point along x in real space. 
However, it remains to determine the dependence of these 
perturbation quantities along the z axis. Thus, eliminating w(k,z) 
and p(k,z) from this set leads to a fourth order ordinary 
differential equation (O.D.E.) with constant coefficients. 



9“*u (k , z) 
9z'‘ 



9^u (k, z) 

(k^ + ikUo/K) 

9z^ 



+ 



ikUpk-u (k , z) 



K 



0 . 



(8.1.9) 

As is the usual fashion for ordinary differential equations, we can 
assume exponential solutions in z of the form, exp(-az), where a is 
a complex number. This leads to the fourth order polynomial 
equation , 
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(k“ + ikUo/K)o“ + ikUok^/K = 0 



( 8 . 1 . 10 ) 



a 



4 



Which is characteristic of the fourth order O.D.E. This equation 
has the four roots, 



a 



+/- jk] , and 
+ /- V(ikUo/K) 



( 8 . 1 . 11 ) 



However, the two positive roots lead to unlimited increases in u 
with z. Thus, only the two negative roots represent components of 
the real solution, 

u(k,z) = U| exp (■/( ikU()/K) z) + U, exp(-|k|z) , (8.1.12a) 

or 

u(k,z) = U| exp(-a, z) + U 2 exp(-a 2 z) . (8.1.12b) 



Note for each component solution that 



3u(k,z)/3z = Gxp(-aj z) 



(8.1.13) 



So in general a is equivalent to the differential operator. 



a = d/dz 



(8. 1.14) 



We will use this below in determining U, and U 2 . But before doing 
so, we focus on the exponential terms. First, we introduce the 
inner and outer length scales, £, and L. In the simplest case, L 
is the horizontal dimension for a solitary hill. This will 

introduce only one wave number, k = 1/L, into the flow response. 
We define the inner scale, €, by 

^ In (f/Zo) = k^L = k , (8.1.15) 

where Zq is the roughness length (Jackson and Hunt, 1979). 
Physically, we are assuming that the hill interacts strongly with 
the region below f through turbulent exchange. Thus, below £ is a 
region of large wind shear, whereas above t turbulent exchange is 
small and resulting in nearly inviscid potential flow. 
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For neutral flow Monin-Obukhov similarity theory specifies for 
horizontally homogeneous surface layers that 



Uq(z) = u,ln{z / Zq) I K , and (8.1.16a) 

K(z) = \i,z/k , (8.1.16b) 

where k, the von Karman constant is = 0.4. Thus, Uq/K = 

In ( z/ Z q) / This leaves, 

a, =V(ikUo/K) = ( 1+1 )a/( In ( z/Zq) / k) A^2 , (8.1.17) 

where we have used a/ 1 = (l + i)/v^2. Now at the height z = £, we have 
cr, = (1 + i) /\/(2£) , (8.1.18) 



while the second root is simply. 



«2 = "1^1 = "1/L • (8.1.19) 

Thus, we have framed the advection velocity component at 
wavenumber, k, in terms of the two Jackson-Hunt length scales. 

Returning to the U, and Uj coefficients, for real flows with a 
non-slip bottom boundary, u(k,z) must become zero at z = 0. From 
eqn. (8.1.12b) this forces 

U, = -Uj . (8.1.20) 

For the magnitude of Uj , we assume to first order that the vertical 
wind near the surface is given by the vertical component of the 
wind vector parallel to the surface, whose speed is taken to be 
equal to the horizontal wind speed. That is, U • Vh = w(0), where 
h is terrain height and U is the steady state wind. By the Fourier 
transformation, h(x) becomes -ikh(k) , where h(k) is defined by 



1 00 

h(x) = J hk exp(-ikx)9k . (8.1.21) 

277 -00 
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This leads to the condition, 



-ikh(k) 



w, 



u 



01 






u, 



02 



( 8 . 1 . 22 ) 



Here we have approximated the mean wind Uq from eqn. (8.1.16a) in 
terms of two components, Uq, and Uqt, as functions of a, and ot 2 ^nd, 
hence, £ and L. Namely, Uq, and Uq, are the mean wind velocities at 
z = £ and z = L. We approximate K, the turbulent diffusivity, 
similarly in terms of £ and L. This implies that a perturbation 
penetrating to larger heights will be subject to larger mean 
velocities and dif f usivities than less penetrative perturbations. 

At any rate, from eqns. (8.1.8) and (8.1.14), we have 



w(k,z) = ik u(k,z)/a . (8.1.23) 



Combining eqn. (8.1.20) with (8.1.22) and (8.1.23) gives 
ik ik 

+ U, = -ikh(k) . (8.1.24) 

°' l^01 ® 2^02 

With LUq, >> £Uq| we can neglect the first term and with eqn. 
(8.1.19)^ 

U, K -lk|UQ2 h(k) . (8.1.25) 



Substituting this into eqn. (8.1.12), we finally obtain for the 
simple 2-D case the velocity perturbation for a single wave number, 
k, at height z above the terrain. 



h(k) 

u(k,z) = Uq(L) [e'^-''" - , (8.1.26) 

L 
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We then back Fourier transform to real space and add the 
perturbation to the mean velocity to produce the total U field. A 
plot of the vertical profile shows that the maximum perturbation in 
2-D occurs at (3 /2v2)£ ' 3.3Z. That is, hills induce pressure 

gradients which interact with the shear stress to force advection 
to produce a low-level jet, as in Jackson-Hunt theory. 

Generalizing to 3-D flow requires that we re-introduce the Coriolis 
term and add an equation and terms for the y component. In the 2-D 
case we solved the O.D.E. system by eliminating variables and 
consolidating equations until only one fourth order O.D.E. 
remained. However, for 3-D cases this procedure is too cumbersome. 
So we re-organize the mathematics by writing the governing 
equations in matrix form and solve them systematically through 
elementary row operations on the matrix, rather than by 
elimination . 

For the 3-D case the form of the Fourier transform and back 
transform now become. 



u (k, m, z ) 



00 00 

f J u (x , y , z) exp ( ikx + imy) 3x dy 

-00 -00 



(8.1.27a) 



1 00 

u(x,y, z) = — J 

277 -00 



00 

j u (k , m, z) exp ( -ikx-imy) ^k 3m 

—00 



(8.1.27b) 



or more precisely, using the discrete (fast) Fourier transform 
pair. 



n-1 n-1 

u(k,m,z) = E E u(x,y,z)W„^^ , (8.1.28a) 

k=0 m=0 



1 n-1 n-1 

u(x,y,z) = E E u(k,m, z)W„.^^ W^_^y , (8.1.28b) 

n^ x=0 y=0 



where W = exp(i27r/n), and similarly for v, w, p, and h. Then the 
momentum and continuity equations for the Fourier transformed 
components at each wave number pair, (k,m), may be written in the 
form of the matrix equation. 



60 



C -f 0 ik 




u (k,m) 




0 


f C 0 im 




v(k,m) 




0 


0 0 C -ct 




w (k , m) 




0 


ik im -a 0 




P(k,m) 




0 



where, using eqn. (8.1.14), 
C = i(kUo + mVo) - Ka“ 



(8. 1.29) 



(8. 1.30) 



Rather than by elimination, we can solve this set systematically 
through elementary row operations (cf. any text on O.D.E. systems) . 
In place of eqn. (8.1.10) the characteristic equation becomes, 

C(k“ + o“) = 0 . (8.1.31) 



We replace eqn. (8.1.15) for f by 



C In (1/Zq) = x:-L/cos (3 



(8.1.32) 



where the mean wind vector, V s (Uq,Vq) , L = 1/k, and B is the slope 
angle, such that 

I V| I (k,m) I cos B = V • (k,m) . (8.1.33) 

The surface boundary condition for w changes from eqn. (8.1.22) to 

W, W2 

-ikh^^ = + (8.1.34) 

(k,m)- Vq, (k,m) • 



For winds not perpendicular to a 2-D hill, Troen (1986) showed that 
eqn. (8.1.26) still gives the perturbation of the component normal 
to the hill. That is, the first, outer scale, term in L is not 
affected by wind direction, but f, given by eqn. (8.1.32), is now 
larger for off-normal winds; below f, the wind speed perturbation 
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is no longer given by a simple cosine response. 

An idealized neutrally stable boundary layer should have no abrupt 
top. In reality the top, H, is usually specified by other 
conditions, such as a subsidence inversion (chronic for Vandenberg) 
or the advected "fossil" remains of a capping inversion formed when 
the upstream boundary layer was last convectively unstable. So we 
must account for H. We also expect that LINCOM will be used in 
cases where the atmosphere is not truly neutral. That is, the 
inverse Obukhov length, 1/L^o will not be zero. In that case, 
can be measured from surface data and used in similarity functions 
in order to obtain the vertical profiles of mean wind speed and 
diffusivity. E.g. , the mean wind and diffusivity profiles may be 
obtained from the following extended boundary layer similarity 
profiles of Sorbjan (1989). He suggested that 



Vq = 2 . 5u.(ln(z/Zo) + 5 z/L - z/H - 2.5 z'/HL^^) ) , (8.1.35a) 



K = 0.4 u.z(l - Z/H)/(1 + 5 Z/L^J 



(8.1.35b) 



for stable cases and 



Vq = 2.5u,(ln(z/Zo) + 0^ 



(8.1.36a) 



K = 0.4 u.z(l - 28 



(8.1.36b) 



for unstable cases. Here we supply a new, more efficient algorithm 
0^ = (1.19 + 0.23*ln(-z/L^J )2 . (8.1.37) 



8.1.4 Obukhov length 

Unless sensors are frequently calibrated, potential temperature 
differences between vertical levels cannot be measured accurately 
enough on a routine operational basis to maintain confidence in 
atmospheric surface layer stability estimates. Thus, LINCOM 1.1 
does not determine Obukhov length from the potential temperature 
and wind speeds differences between vertical levels. Instead, it 
relies on measured solar and thermal radiation values or estimates 
thereof. This type of method is useful for Vandenberg because the 
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near coastal areas are chronically shrouded with stratus cover, 
while the sky is often quite clear several kilometers inland. Thus, 
for diagnostic models such as LINCOM, Obukhov lengths need be 
computed for only a few areas between sun and shade, based on 
moderate differences in roughness length and local wind speeds. 

We also included a routine to estimate downwelling solar and 
thermal radiation, if measured data is not available. This routine 
requires a fractional cloud coverage estimate as well as other 
commonly available input data, such as inversion height, screen 
level wind speed, temperature, and relative humidity. 

All solar radiation models require current date and time to compute 
sun angle, \p. This is done in LINCOM using standard algorithms. We 
modified a simple ASHRAE model for downwelling solar direct and 
diffuse radiation (Iqbal, 1983) by computing sine curve fits to the 
seasonal adjustments for apparent extra-terrestrial radiation, 
atmospheric optical air mass, and column length of precipitable 
water, W. 

If cloud base height, and boundary layer height, H, are not 
available, the solar transmissivity through clouds, is 
estimated using the simple algorithm. 



= (1 - 0.6FCC) 



(8. 1.38) 



where FCC is fractional cloud coverage of the sky. FCC can be 
obtained from ground based observation and/or on-line GOES 
satellite data, using the algorithm described in Skupniewicz et al, 
(1991). If Zj, and H are available (from rawinsonde data or 
aircraft landing reports), the solar transmissivity for a non- 
reflective ground surface can be estimated more accurately from the 
method of Liou and Wittman (1979). This method uses sun angle and 
column height of precipitable water within the cloud in a bivariate 
polynomial regression based on results from an accurate multi- 
stream discrete ordinates model. Currently, we assume that the 
cloud coverage is stratiform and confined to the boundary layer, 
with a liquid water content of 0.78 grams/meter of cloud depth. 
Hence, the only inputs required are date, time, cloud base height, 
and boundary layer depth. The algorithm can be extended easily to 
include other cloud types and water content. The regression form is 



Tc(Mo,W) 



3 3 

S bjj M'o W' 

i = 0 j = 0 ■' 



(8.1.39) 



where is solar zenith angle, W is precipitable water, and the b- 
are the coefficients obtained from the regression. For actual non- 
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zero surface albedoes, (default value, 0.15) , we modify the cloud 
solar transmissivity, using the algorithm of Kamada (1984), 



FCC 



1 - A,A, 



^cm 



(1 - 0.12A^) (1 - O.lAj) 



(1 - d0)T. + d<^) 



(8.1.40) 

where A^. is cloud top albedo ( default value for stratiform clouds 
is 0.55), d = 0.001068, and <p is in degrees latitude. The total 
downwelling solar radiation is then 

SOLI = (IoSin(i/^) + I,,iff)T,^ , (8.1.41) 

where is the downward diffuse solar radiation at the surface. 

The downwelling thermal radiation computation is initiated, using 
the algorithm of Martin and Berdahl (1983) , which employs surface 
dewpoint temperature, T^j^, hour of the day. Hr, and pressure, p, to 
estimate the effective clear sky emissivity, 

€„ = 0.711 + 0.0056 Tjp + 0.00073 T^p^ + 

0.013cos(0.262Hr) + 0.00012(p - 1000) . (8.1.42) 



Tjp is readily obtained from the relative humidity and temperature, 
using standard formulas. The emissivity is then modified for clouds 
according to cloud base height, z^, and fractional cloud coverage, 
FCC. Thus, we have 



= €„ + 0.85FCC(1 - e„) exp (1. 22x10'' Z^) . (8.1.43) 



Again, boundary layer stratus clouds are assumed here but other 
cloud types are readily included. Downwelling thermal radiation is 
computed by assuming the cloudy or clear sky to be a grey body 
thermal emitter, such that 



IRl = €^.09'' 



(8.1.44) 
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where a = 5.67x10® is the Planck black body constant. Total 

downwelling radiation at the earth's surface is obtained by 
combining solar and thermal contributions via, 

RADI = (1 - Aj,)S0Ll + IRl . (8.1.45) 



With the radiation component of the surface energy budget computed, 
we can obtain the atmospheric stability and temperature flux from 
the ground surface to the air. As a measure of atmospheric 
stability, the Obukhov length is defined as, 

L = u.-e/gxG. , (8.1.46) 

where u* is the surface layer friction velocity, 6 = T ( 1000/p) 
is the near surface potential temperature, g is gravitational 
acceleration, k is 0.4, the von Karman constant, and 6» is the 
Obukhov temperature scale. The temperature flux can be defined as 
the statistical correlation between vertical velocity and potential 
temperature perturbations, and is also given by 



w '0 *Q = -u. 0. 



(8 . 1 .47) 



Thus, given the downwelling radiation, we can iterate between 
estimates of L and estimates of surface temperature flux until both 
quantities converge. This requires that we compute both u» and 0*. 
u. comes from 



u, = U(z) k 

J Ln(z/ZQ) - '^m 



(8.1.48) 



where U(z) is the mean windspeed at height z, and Zq is the surface 
vegetative canopy roughness length (typically = 1/7 the mean 
vegetation height). 0^ is given by eqn. (8.1.37). Analogous to u,, 
0, is given by 



0, = <50(z) I k 

] Ln(z/Zo) - \ 



(8.1.49) 
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From Dyer and Bradley (1982), we have 



= 2 Ln ( ( 1 - \ 1 - 142/L ) ) 



(8 . 1.50) 



To obtain the ground surface "skin" temperature, ©qq/ we assume that 
the potential temperature difference between the roughness height, 
Zq, and height, z, is given by 60, and that the temperature 
difference across the laminar layer between Zq and the surface is 
given by ©,. Thus, 



©00 = e - (©./O (Ln(z/Zo) - \) - ©. . (8.1.51) 



Since 6 © is not known, ©. is initially set to zero and the skin 
temperature is set equal to the screen level potential temperature. 
The temperatures then diverge toward equilibrium values by 
iteration. The net radiative budget at the surface is given by 



NETRAD = RADI - Cg o ©oo"* , (8.1.52) 

where the ground surface emissivity, € , has a default value of 
0.95. Following the Penman-Monteith model (1948, 1965), 
temperature flux into the ground is estimated as, Q = 0.15 NETRAD 
with the stipulation that for stable conditions (if L > 0) , Qg is 
3.3 times larger. The Penman-Monteith equation is used to estimate 
the temperature flux. 



-(7 (-NETRAD + Qg) -w'q'o) 

w'©'o = - 0.84w'q'o, 

( Q Cp (XgS^^ + 7) ) 



(8.1.53) 



where w'q'p is the humidity flux, p is air density, = 1005 J kg'* 
K'* is the heat capacity of air, is soil relative humidity, s^^ is 
the change rate of specific humidity with temperature for saturated 
air, and 7 = Cp/L^, = 0.0004°K* is the psychrometric constant. In 
turn these latter variables are obtained from standard algorithms. 
The Obukhov temperature scale is obtained from. 



©» = -w'©'q/u» 



(8.1.54) 
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In this scheme, note that under neutral conditions when the 
temperature flux falls to zero, 0« will vanish and L becomes 
infinite. Thus, to avoid infinities during iterative numerical 
evaluation, it is better to compute the inverse Obukhov length, 
1/L. 

From a rough bivariate analysis of our results, we found that a 
useful first guess is 



1/L = 0.000674 (300 - RADi) Ln ( 10 z / Zq) /U‘ ( z ) . (8.1.55) 



In order to cover a wide range of radiation values, wind speeds, 
temperatures, and roughness lengths, the iteration procedure must 
be quite robust, otherwise convergence is not obtained, especially 
at low wind speeds. Even with a good first guess like eqn. 
(8.1.54), we found that simple iteration was unreliable, that the 
Newton-secant root finding procedure was needed for standard cases, 
and that second order Aitken acceleration was required for low wind 
speeds and small roughness lengths. 

The Newton-secant root finding algorithm utilizes the form. 



+ , - Xj = <Si^, = 



-<5, f (X;) 
f(Xj) - f(Xi_,) 



(8.1.56) 



where i is iteration number. Here x is 1/L, and f(x) is the 
difference between old and new values of 1/L. The object of the 
iteration process is to adjust f(x) toward zero. Unlike the 
standard Newton-Raphson technique, the secant method does not 
require an analytic expression for the first derivative of f(x), 
which in our case is not obtainable. The Aitken technique is a 
second order acceleration method. 



XiXi + 2 - X i+i 

Xj + 3 “ Xj = = , 

Xj + 2 ~ 2Xj+, + Xj (8.1.57) 



which relies on the Newton-secant method for the first two 
iterations, then computes the rate of change of the convergence 
from the first two iterations and uses it to obtain a refined 
estimate for the third and subsequent iterations. 
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8.1.5 Objective Analysis 



To produce each of the n flow fields for the n towers, the mean 
wind vector from each tower is added to the terrain induced 
perturbation field from LINCOM. Then at each grid point, the Uj 
from the n flow fields are combined according to the algorithm, 

n 

Z Ui/r;^ 

i=l 

U(x,y) = , (8.1.58) 

n 

Z l/r;* 
i=l 

where rj is the terrain modified distance from the grid point to 
tower i. In general we define r as. 



r = (a + (b - a)|sin 0|)R 



(8.1.59) 



Here <p is the tower-to-gr id-point angle measured from north (as in 
meteorological wind angle) . R is the actual tower-to-grid-point 
distance , 



a = V{h/B) , and b = 1/a 



where the dominant aspect ratio of the surrounding terrain height 
contours is given by A/B. Thus, rather than a 1/R^ diminished 
circle of influence, each tower has a 1/r^ diminished oval of 
influence which corresponds to the dominant shape and orientation 
of the elevation contours surrounding the tower. Of course for 
every grid point the rj in a given domain need be computed only 
once, then stored for future use. 

This method of combining flow fields from each tower ensures that 
the final field will exactly match the tower data, while winds at 
locations between towers are prescribed by a combination of terrain 
modified objective analysis and neutral boundary layer dynamics. 

In LINCOM 1.1 we noted from eqns . (8.1.16 - 19) that aj and aj, 

the vertical decay coefficients for the terrain induced wind 
perturbation fields, are independent of mean wind speed. For 
neutral flow this means that the fundamental characteristics of the 
perturbation field are determined by the terrain and are thus 
domain specific. So the perturbation field may be expressed as a 
linear combination of the mean wind vector and the two orthogonal 
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component solutions for the perturbations along the x and y axes. 
Since the component solutions for any given domain may be 
pre-stored, this allows an extremely fast "look-up table" mode of 
running LINCOM. 

Actually, several sets of orthogonal solutions may be stored to 
account for different inversion heights. Interpolation and a 
fitting algorithm are then used to determine the mean wind vector, 
which when combined with the perturbation fields, will best fit the 
available tower data. The best fit is determined by minimizing the 
differences between the tower and LINCOM winds. Unlike LINCOM 1.0 
this new procedure in LINCOM 1.1 relaxes absolute tower anchoring 
in order to ensure mass conservation of the flow field. A 
reasonably close fit to the tower data is retained, but with little 
potential for aliasing as described in sections 3.1 and 5.2. This 
procedure also results in a significant gain in computational speed 
when many towers are involved. 
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8.2 RIMPUFF Model Theory 



8.2.1 Introduction 

Since distinctions between turbulence and mean wind are not well 
resolved in complex terrain, some arbitrary scale such as the 
model's grid resolution becomes a convenient divider. That is, we 
can determine the mean flow down to the grid scale as a direct 
deterministic response to the physical forces. However, the grid 
resolution cannot be arbitrarily fine because we lack fine scale 
information about bottom boundary conditions, especially in complex 
terrain, and also because computer power is finite. Moreover, 
sub-grid scale effects must be parameterized rather than simply 
neglected because the non-linear eddy stress dynamically links all 
scales of turbulence in an energy cascade. Thus, without modeling 
the dissipation to heat, the turbulence kinetic energy would simply 
grow without bound. Luckily for our investigation, the sub-grid 
fluid motions can be treated statistically, since smaller scales 
tend to behave more randomly, even in complex terrain. 

Hence, RIMPUFF lets LINCOM handle the bulk transport, while it 
treats the instantaneous diffusion in a relative frame which 
follows the mean flow. Continuous releases, i.e., plumes, are 
modeled by a series of puffs released into the meandering flow 
field. RIMPUFF computes species concentrations by summing 
contributions from all puffs for each grid point and time step. 
Depending on computer capacity, the model can handle several 
hundred puffs released from multiple point sources anywhere within 
the domain. Such sources can have individual release rates, times 
and heat production rates. Plume heights depend on standard 
formulas using downwind distance, stability, and stability 
dependent, vertical profiles of windspeed. The inversion is 
assumed to parallel the terrain and is treated as a material 
surface. Both the inversion and source heights are user specified. 

Total surface reflection of puffs is normally assumed for inert 
gases, but the reflection coefficient is adjustable from zero to 
unity for materials with real deposition rates. Dry deposition is 
calculated from source depletion. Dispersion parameters for each 
puff depend on items such as pollutant type, atmospheric stability, 
and wind speed. Wet deposition depends on species and the space/ 
time distribution of the rain field. 



8.2.2 Governing equations 

The horizontal sub-grid diffusion is treated by linking puff growth 
to the time averaged, local component spectra of the lateral 
velocities (see below) . A kinematic-statistical approach is 
employed which assumes that the relative displacement of fluid 
elements proceeds in a gaussian manner. 
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That is, with Q referring to total puff mass and C(x,t) to the 
observed concentration field in space and time, we begin by posing 
the mass conservation and center of mass location as 



Q= Jc(x,t)5x , and 



( 8 . 2 . 1 ) 



c(t) = 1/Q Jx C(x,t)3x 



( 8 . 2 . 2 ) 



The center of mass is subject to transport by the mean flow 
determined by LINCOM and also by grid scale and sub-grid scale 
turbulence not included in LINCOM. This turbulence "meander” is 
estimated by a stochastic puff advection scheme described later. 
At any rate the center of mass velocity in a fixed frame can be 
described as a summation of the product of the puff mass at each 
location with its respective velocity. 



where u(x,t) represents the fluid velocity field. We then denote 
both the average of all particles within a puff and the ensemble 
average by "< >" and fluctuations by " ' ". Thus, if C(x,t) = 

< C(x,t) > + C’(x,t), then 1/Q< C(x,t) >dx denotes the probability 
of finding a marked particle at the point (x,t) and 

< X^ > = 1/Q J x*< C(x,t) >dx (8.2.4) 

denotes the ensemble averaged, fixed frame, absolute diffusion. 
The latter may be expressed as the sum of meander and relative 
diffusion, or 



< x2 (t) > = < y^ (t) > + < c^ (t) > , (8.2.5) 

where y = x - c defines a coordinate frame moving with the center 
of mass velocity. 

Now, using t to time integrate along a particle’s path with 
respect to the center of mass, the puff growth rate becomes, 



V^^(t) = 1/Q J u(x,t) C(x,t) dx 



(8.2.3) 



t 

da^/dt = 2\ < v(t) v(t-T) > dT 



o 
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< (t) > t,(t) 



t 

= < (t) >J r(t,T)dr = 

o 

( 8 . 2 . 6 ) 

Here, a is puff size (standard deviation of particle displaceroent 
relative to center of mass) and < v^ (t) > is the corresponding 

relative velocity variance. The relative correlation function, 
r(t,T), depends on both t and T. This indicates a non-stationary 
process, unlike Taylor's single particle diffusion. Letting v = u 
- be the moving frame particle velocity, we can show that 



t 

do^/dt = 2 J [Rahs(0 - Rcn,(t,r)) dr 
o 



(8.2.7) 



f{ere , 



Rabs(’’) = < u(t)u(t-T) > , and 



(8.2.8a) 



Rcm(t,r) = < V,^(t)V,^(t-T)) > 



(8.2.8b) 



are the appropriate auto-covariance functions for absolute single 

particle diffusion, and from the center-of-mass-velocity viewpoint, 

respectively. Since Raj,,, is known, we can focus on obtaining a model 

for R-_. 
cm 

Using (8.2.3), R^.^(t,T) can be rewritten as 



n < u (x ' , t) u (x" , t-T ) C (X ' , t ) C (x" , t-T ) >3x'3x" . (8.2.9) 

From this point the puff can be modeled as a set of passive marked 
particles with instantaneous concentration, C(x,t), moving within 
a large but numerable set of unmarked fluid elements, with Eulerian 
velocity field u(x,t). Any relative displacements, y(t), can also 
be assumed to obey the same un-correlated gaussian statistics (see 
fig. 8.2.1). That is, any two particles viewed in the moving frame 
will appear to disperse in an uncorrelated fashion, such that 



<Ay Ay > = 

i j 



(t) - a* (t-T) for i = j 

0 for i j 



( 8 . 2 . 10 ) 
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Physically, this means that the characteristic scale of motion in 
the moving frame must be much shorter than the puff size. However, 
these scales may overlap in the atmosphere, implying that only the 
ensemble average of individual puff releases will be approximately 
gaussian . 

Hence, a gaussian puff of initial size, a(0), will remain a 
gaussian puff of size, a(tj), at a later time, (tj) , so that 

= J G5(y-yo,t) G„(,.,)(yo,t-T) 3yo . (8.2.11) 



We also see that the spread, < Ay*j(t,T) >, for the transition 
probability, Gj, given in eqn. (8.2.10) satisfies eqn. (8.2.11). 

Substituting the instantaneous concentrations from eqn. (8.2.9) 
with the corresponding gaussians now leads to 



Rcm(t,0 = n< u(y' ,t)u(y",t-T) >G^(,)(y’ ,t)G„(,,,,(y",t-T)ay'ay”, 



( 8 . 2 . 12 ) 

where the moving frame fluid velocity, u, relates to the fixed 
frame velocity by u(y,t) = u(y+c,t). 

We must also obtain the two-point, two-time (fixed frame) velocity 
covariance, < u (y ' , t) u (y" , t-r ) >. For two particles, i and j, in 
a homogeneously turbulent field, the velocity covariance function, 
in a fixed frame, (^jj , t ) , will depend only upon the temporal and 
spatial separations, t and between the particles. That is. 



= < u(x,(t)) u(Xj(t-T) - ^ij) > . (8.2.13) 

Thus, the moving frame velocity covariance in eqn. (8.2.12) is 
written in terms of the absolute velocity covariance function and 
the transition probability from eqn. (8.2.11) as 



< u(y’ , t)u(y",t-T) > = i Gj(y'-y"-C,t)Raj„(^,T)d^ . (8.2.14) 

Returning this to eqn. (8.2.12) and integrating over y* and y" will 
reformulate the center-of -mass covariance as 
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(8.2.15) 



Thus, in terms of the fixed frame covariances, the horizontal 
rate equation for puff growth becomes. 



t 

aavat = 2 \ (Rabs(o^o - 

o 

RabJ^O/t 2 v/tt a(t) d^ar 

(8.2.16) 



For ease of use, we employ the Fourier transform. 



R,b,(^T) = S(«,w) e‘('<^+“^)a^aw 



(8.2. 17) 



to render eqn. (8.2.16) as. 



da^/dt = 2tJJ S(»c,cj) sin(ojt)/ojt (1 - dc^dK 



(8.2.18) 

In this case the spatial filter, 1 - has been inserted to 
remove wave numbers smaller than 1/a, while the low pass temporal 
frequency filter, sin (wt) /cjt , does likewise beyond o) = 1/t. The 
shaded area in Fig. 8.2.2 shows the part of S(k,cj) that effectively 
contributes to the growth rate. 

For puffs much larger than the turbulence length scale, eqn. 

(8.2.18) reduces to Taylor's formula for single particle diffusion. 
Hence, behavioral differences between a puff and a single particle 
are closely related to the spatial correlation of the turbulence. 
That is, for classically homogeneous turbulence with a Kolmogorov 
inertial subrange, the spectral index, p = -5/3. Puff growth will 
then follow a standard prediction. However, in a stably 
stratified atmosphere, p may be = -3, implying exponential growth. 
In either case, we assume that the sub-grid diffusion is locally 
homogeneous. If so, we can scale the puff's instantaneous growth 
rate over intermediate ranges (< 500m) with the local lateral 
turbulence intensity, a;, such that da/dt « 0.3ajU. We obtained 
turbulence intensities for Vandenberg by modeling Eulerian wave 
lengths with the similarity formula. 



0; = ajo( 1 + B/ut)-'^3 



(8.2. 19) 



where refers to limiting values at large times ( Skupniewicz , et 
al, 1989). That is, even within complex terrain, the turbulent 
energy characterized by Vandenberg's tower wind spectra tends to 
plateau at larger scales. Equation (8.2.19) tends to apply over a 
wider range of scales than simple power law estimates. Power law 
fits to Vandenberg data are also found to be rather windspeed 
dependent. Note that the constant, B, depends upon measurement 
height and stability and is proportional to the dominant eddy size 
at a given height. Thus, values for Cjq and B were determined by 
multi-variate regression for each tower for each of three periods: 
pre-dawn, noon, and dusk, indicative of stable, unstable, and 
neutral conditions. Values were also classed by wind direction. 
The longitudinal and lateral values determined for o-^q ranged 
typically from 700 - 2,000m and 200 - 1,000m, respectively. As in 
flat terrain, the longitudinal component clearly dominates. 
However, our values were consistently larger, suggesting that flow 
distortions due to complex terrain feed energy directly into larger 
horizontal scales of motion. 



8.2.3 Puff splitting 

Real diffusion in complex terrain also often displays plume 
bifurcation and/or vertical shear due to channeling, slope flows, 
or inversions. To simulate such decoupling we allow each puff to 
split as often as necessary. The daughter puffs are arranged such 
that the original center concentration, radially integrated second 
moments of the concentration distribution, and mass are all 
conserved . 

These constraints may be summarized as. 



2 



i 




( 8 . 2 . 20 a) 



5 



ii € 5 ( 0 , 0 ) = C,(0,0) 



( 8 . 2 . 20 b) 






(l/2)ap, 



and 



( 8 . 2 . 20 c) 



IV 



mass conservation. 
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For the Vandenberg simulations the initial puff diameter was set to 
100m ( one standard deviation in concentration) . When the puff 
attains a diameter equal to the 500m grid spacing, it splits 
(pentafurcates) horizontally into five new gaussian puffs. The 
daughter puff diameters are initially set to 250m. The 
conservation relations then require that the four satellite puffs 
be centered at 0.89 from the origin and that the satellite and 
center puffs each carry 23.53 and 5.88% of the total matter, 
respectively (fig. 5). 

In the case of layers decoupled vertically by shear, the original 
puff is allowed to trifurcate into three daughter puffs which are 
centered along the original puff's vertical axis. To simulate a 
gaussian radial concentration profile, the conservation rules 
require that the two satellite puffs be centered at +/-1.19a^, with 
each new puff carrying 26.58 % of the initial pollutant mass. 

The new center puff retains the remaining 48.84 %. Again, the new 
puffs are born with half the diameter of the original puff. 

The purpose of this splitting scheme is to allow each new puff its 
own individual growth rate and direction, as set by the mean flow 
and turbulence intensities local to it. We assume in this 
procedure that the sub-grid scale shearing is locally homogeneous 
so that filtered turbulence estimates can be used in the 
parameterization. The puff-splitting scheme can be repeated to a 
practical limit of several hundred puff progeny. In this way the 
plume bifurcation caused by terrain features can be modeled in 
detail. Such bifurcation was seen extensively during LVDE (see 
LVDE Data Report, 1990) 



8.2.4 Stochastic puf f-advection 

The mean flow field provided by LINCOM should be updated every 5 to 
10 minutes to include temporal variability and allow for puff 
meander in RIMPUFF due to non-stationary winds. But if measured 
winds are available only as one hour means, only horizontal shear 
in the static wind field is left to induce any meander in the puff 
center-of-mass advection velocity, V,,^. Such was the case for the 
eight Mt. Iron diffusion experiments. To re-introduce temporal 
variability in the wind field, we have added an auto-regressive 
(Langevin equation based) advection component to the two horizontal 
(u, v) wind components of the mean flow model. This is a)cin to 
standard Langevin random forcing in Monte Carlo dispersion but 
applied to puffs rather than particles. 

We estimated supra-grid scale (> 500 m) wind variance, by 
low-pass filtering the horizontal wind energy spectrum below wave 
number k = 2a/500 m. Since we already use the sub-grid scale (< 
500 m) variance for puff growth, we now account for turbulence at 
all scales. Tl p„ff = L^/u also provides a Lagrangian time scale for 
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horizontal puff motion, where Lg, an Eulerian length scale, is 
derived from velocity spectra measurements (see Handbook) . 

With these we can estimate stochastic contributions to from 






( 8 . 2 . 21 ) 



where the Lagrangian correlation coefficient for the center-of-mass 
velocity is given by 



p = exp(-At/TL) , (8.2.22) 

and rj is white gaussian noise having a variance, (1 - P“) . 

By comparing cases with and without, we found that the stochastic 
component is clearly visible during individual simulations. 
However, its contribution to the overall plume widths is small. 
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TIME 



Fig. 8.2.1 The moving frame trajectory, y;, of a 
particle, i, that at time t-r , holds the position, 
guantity, Ay; = yj(t) - yj(t-T), and its gaussian 
function, G^, are shown at time t. 



marked fluid 
yi(t-T). The 
distribution 
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Fig. 8.2.2 Schematic isopleth plot of spectrum, S(k,w). Its 
maximum value is at = (0,0) from where the function decreases 
monotonically through the levels I, II, and III. The shaded region 
is the effective domain of the model defined by vertical line, h = 
f*, (using the low-pass filter, sin(bt)/bt) and the horizontal 
line, k = a* (using the high pass filter, 1 - exp(-k^a^) ). 
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Peixia Puff 




Fig. 8.2.3 Puff pentaf urcation : the original single puff of size 
a divides into five new puffs, each of size h t conserving total 
mass, and peak and second moments of the concentration distribution 
around the cloud's center of mass point (moment of inertia) . 
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8.3 Comparison Procedures 

The following describes how the LINCOM/RIMPUFF model runs were 
compared with the digitized Mt. Iron data for cases 28, 31, 48, 55, 
87, 90, 91 and 110. 

First, the LINCOM/RIMPUFF model runs and Mt . Iron case studies were 
plotted with the SURFER graphics package so the plumes could be 
visually inspected. By overlaying a model run with a case study it 
was evident for some cases that the source points were not 
collocated. The Mt. Iron grids were moved by an appropriate number 
of grid points in the x and y directions to align the source point 
of the Mt. Iron case study and the LINCOM/RIMPUFF model run. 

Headers at the top of the SURFER formatted *.grd files define the 
latitude/ longitude coordinates in UTMS for the four corners of the 
LINCOM/RIMPUFF and Mt . Iron dose files. For all of the LINCOM/ 
RIMPUFF runs, the headers did not correctly locate the source 
points relative to a map of Vandenberg. By plotting some of the 
roads at Vandenberg along with the location of VIP-1, the proper 
coordinates for the surfer header files for each of the model runs 
were determined. These roads were plotted by picking off latitude 
and longitude coordinates for points on the roads using a 300' (.3 

arc seconds) resolution map of Vandenberg. Program utm.for was 
written to convert the longitude/ latitude values to UTMS 
coordinates . 



For cases 28, 48, 87, 90 and 110 the plotted roads (which were also 
sampling lines for the Mt. Iron data set) were used to define plume 
cutoff lines. Since the plumes did not reach Santa Ynez Rd . in 
cases 28 and 110, these plumes were truncated at Arguello Rd . The 
plumes were truncated at Santa Ynez Rd. for cases 87 and 90 because 
there was enough data to define these plumes at Santa Ynez. For 
case 48, the plume was cut off at Honda Ridge. 



The following procedure was used for this evaluation. 1) TEST, an 
ASCII file was created to hold the input and output names of the 
LINCOM/RIMPUFF and Mt . Iron grid files. This file is comprised of 
four file names followed by a blank line for each comparison. For 
example, to compare a digitized data grid with a model simulation 
grid for case 28, TEST would look like: 



mi28 .grd 
case28 . grd 
p28 .grd 
c28 . grd 



# input file name of Mt. Iron data grid 

# input file name of LINCOM/RIMPUFF simulation grid 

# output file name of Mt. Iron grid 

§ output file name of LINCOM/RIMPUFF grid 



2) The program, drawsurfl.f, was run to collocate the release 
points, find the merit scores where the hand-drawn isopleths were 
left untruncated (cases 31, 55 and 91), and to create an output 

file to contour with TOPO . 
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3) The header files in all of the LINCOM model files were changed 
so that the x and y are in UTMS coordinates: 

4) Cases 31, 55 and 91 were then complete and ready to be plotted. 

5) Program mt3.for was run for the truncation cases: 28, 87, 90 and 
110. All concentrations beyond the last sampling line where 
appreciable doses were measured were set to zero. 

6) mt4 . for was run for case 48 which was truncated at Honda Ridge. 

7) Program drawsrf2.f was run on cases 28, 48, 87, 90 and 110. 
This program is actually the same as drawsrfl.f except that the 
format statements to read the input grids are different and also 
the source point adjustment is not re-evaluated. 

8) Cases 28, 48, 87, 90 and 110 were then complete and plotted. 
The "roads. bln" file was used to overlay some roads on the plots. 
Some of these roads correspond to sampling lines. 
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